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Abstract 

\  : 

-l 

The  ultimate  mechanical  properties  of  polymers,  assuming 
perfect  morphology,  will  be  limited  by  the  mechanical  properties  of  a 

r 

single,  ideal  polymer  chain.  Previous  calculations  of  polymer  chain 
moduli  using  semi-empirical  (SE)  quantum  chemistry  methods  have 
resulted  in  modulus  values  much  higher  than  experimentally  measured. 
This  study  investigated  the  error  in  the  calculated  inherent  to  the  method 
of  calculation  by  comparing  SE  results  for  C-C  bond  potentials  in  two  well 
characterized  polymers,  polyethylene  and  polydiacetylene.  It  was  found 
that  the  SE  calculation  systematically  overpredicted  bond  stiffness  in  these 
polymers  by  approximately  25%  to  30%.  This  is  the  upper  limit  on  the 
modulus  overprediction,  depending  on  the  importance  of  bond  exten¬ 
sion/compression  (as  compared  to  other  deformation  modes)  in  the  overall 
deformation  of  the  polymer  chain.  Jt  is  believed  that  this  discrepancy  is 
caused  in  part  by  the  omission  of  bond  stiffness  information  in  the 
parameterization  of  the  SE  method,  and  in  part  by  the  omission  of  electron 
correlation  in  the  calculation  (Hartree-Fock  Limit). 


i 
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I  INTRODUCTION 


The  wave  mechanics  formulation  of  quantum  theory  derived  by  Erwin 
Schrodinger  in  1926  satisfied  the  three  fundamental  requirements  of  a  good 
scientific  theory:  agreement  with  and  explanation  of  existing  experimental 
data,  and  prediction  of  phenomena  which  were  confirmed  in  subsequent 
experiments.  The  power  of  wave  mechanics  is  tempered  by  its 
fundamental  limitation,  that  analytic  solutions  for  wave  functions  and 
energies  can  only  be  obtained  for  atoms  with  symmetric,  central  Coulomb 
potentials,  namely  atomic  hydrogen  (and  one-electron  ions).  One  could 
argue  that  being  able  to  describe  over  99%  of  the  known  matter  in  the 
universe  is  sufficient.  However,  since  the  remaining  1%  includes 
6-electron  atom  based  humans  who  breathe  18-electron  molecules,  there 
was  sufficient  motivation  to  develop  approximate  methods  to  extend 
Schrodinger’s  theory  to  multi-electron  systems. 

We’ll  consider  two  quantum-based  approximations,  both  of  which  rely 
on  the  Hartree-Fock  /  self  consistent  field  (SCF)  method.  The  ab-initio 
implementation  of  the  SCF  makes  few  simplifying  assumptions,  and 
provides  the  best  agreement  with  experiment.  However,  ab-initio  methods 
are  computationally  expensive  and  impractical  for  all  but  fairly  simple 
systems  (roughly  fewer  than  10  non-hydrogen  atoms)  even  on  current  state- 
of-the-art  supercomputers.  Another  implementation  of  the  SCF  method  is 
called  the  semi-empirical  (SE)  approach  because  it  relies  partly  on 
experimental  data  from  atomic  spectra.  The  main  difference  between  ab- 
initio  and  SE  methods  is  that  some  of  the  electron-electron  interaction 


integrals  are  replaced  by  a  set  of  parameters  which  have  been  selected  to 
provide  the  best  agreement  (minimum  least  squares  error)  to  equilibrium 
molecular  properties  of  some  large  set  of  test  compounds.  SE  results  have 
poorer  agreement  with  experiment  than  ab-initio  results,  but  the 
parameterization  decreases  the  computational  cost  so  that  largei  systems 
(up  to  approximately  50  non-hydrogen  atoms)  can  be  readily  calculated. 
However,  SE  methods  are  limited  to  those  elements  which  have  been 
parameterized.  A  third  method,  known  as  molecular  mechanics,  is  based 
on  classical  mechanics  and  has  been  used  on  very  large  systems  (up  to 
100,000  non-hydrogen  atoms).  Application  of  molecular  mechanics  is 
limited  to  situations  where  system-specific  interaction  potentials,  usually 
atom-pair  potentials,  have  been  measured.  Although  this  is  a  non 
quantum-based  method,  molecular  mechanics  provides  excellent 
agreement  with  experiment  because  the  measured  potentials  represent  the 
net  result  of  all  the  physical  processes  involved,  however  small  (i.e.,  makes 
no  approximations).  The  accuracy  of  the  results  depends  directly  on  the 
experimental  accuracy.  These  three  approaches  are  complementary, 
providing  access  to  different  (.not  orthogonal)  parameter  spaces. 

The  semi-empirical  method  has  recently  been  applied  to  calculate 
mechanical  properties  of  polymers  away  from  equilibrium!,2.  In  these 
calculations,  an  infinite,  ideal  polymer  is  simulated  by  calculating  a 
molecule  composed  of  one  or  more  monomer  repeat  units  (a  cluster) 
specified  with  periodic  boundary  conditions.  The  number  of  repeat  units  in 
the  cluster  is  determined  by  the  conjugation  length  in  the  polymer,  i.e.,  the 
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cluster  must  be  long  enough  so  that  its  ends  are  electronically  isolated.  The 
minimum  energy  geometry  of  the  cluster  is  calculated,  with  the  overall 
cluster  length  specified  as  a  variable  which  is  optimized.  This  length  is 
increased  (  or  decreased)  and  then  held  fixed,  and  the  rest  of  the  cluster 
geometry  is  optimized  under  this  constraint.  The  minimum  heat  of 
formation  of  the  cluster  is  calculated  for  a  range  of  cluster  lengths, 
generating  r  restoring  potential  for  the  cluster.  The  second  derivative  of 
this  potential  is  the  effective  force  constant  for  the  cluster,  from  which  the 
polymer  Young's  modulus,  stiffness  per  unit  area,  can  be  calculated.  The 
modulus  values  calculated  using  this  procedure  are  significantly  higher 
than  experimental  measurement.  There  are  two  possible  reasons  for  this 
result:  measurements  are  made  on  bulk  materials,  which  deviate  from  the 
ideal  structure  simulated  in  the  calculations  (ultimate  modulus  vs.  practi¬ 
cal  modulus),  and/or  the  semi-empirical  method  parameters  are  not  valid 
away  from  the  conditions  for  which  they  were  derived. 

The  goal  of  this  thesis  is  to  perform  a  critical  assessment  of  the 
applicability  of  the  semi-empirical  (SE)  quantum  approximation  for  the 
calculation  of  polymer  mechanical  properties  away  from  equilibrium,  and 
to  bound  the  accuracy  of  the  calculation  of  molecular  bond  stiffnesses  at 
equilbrium.  The  non-experimental  parameters  used  in  the  SE  methods  are 
obtained  by  adjuring  a  set  of  initial  estimates  of  the  parameters  for  each 
element  to  minimize  the  least  squares  error  in  reproducing  various 
molecular  properties.  In  the  approximation  used  in  this  study,  the 
parameters  were  optimized  to  reproduce  gas-phase,  equilibrium,  297  K 
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values  of  heats  of  formation,  dipole  moments,  ionization  potentials,  and 
molecular  geometries,  none  of  which  depend  sensitively  on  the  derivatives 
of  the  intra-molecular  potentials.  The  modulus  calculation  involves  a 
series  of  SE  calculations  in  which  one  parameter,  the  length  of  the  cluster 
along  the  polymer  axis,  is  fixed  at  non-equilibrium  values.  All  other 
cluster  properties,  bond  angles,  bond  lengths,  are  adjusted  away  from  their 
equilibrium  values  to  minimize  the  heat  of  formation  under  this  constraint. 
Hence  the  modulus  calculation  involves  extrapolations  of  non-equilibrium 
properties  based  on  equilibrium  data  which  provide  no  guidance  for  the 
direction  this  extrapolation,  i.e.,  derivatives  of  the  intra-molecular 
potentials  .  Therefore,  the  range  of  displacements  away  from  equilibrium 
for  which  SE  calculations  may  be  accurate  is  undefined.  The  equilibrium 
bond  stiffness  is  the  second  derivative  of  the  bond  potential  at  the 
equilibrium  bond  length.  These  data  were  not  used  to  anchor  the  SE 
parameterization,  thus  the  accuracy  of  bond  stiffnesses  predicted  by  a  SE 
calculation  requires  investigation. 

This  assessment  of  the  SE  approach  will  be  performed  by  comparing 
calculated  and  observed  properties  of  polyethylene  (PE)  and  polydiacetylene 
(PDA).  These  polymers  were  selected  because  they  possess  simple  backbone 
structures,  both  have  been  well  characterized  experimentally,  and  they  can 
be  made  in  crystalline  form  which  removes  much  of  the  morphological 
uncertainty  in  comparing  the  experimental  results  to  the  calculations. 

Real  polymers  are  at  best  semi-ordered,  have  a  random  chain  orientation, 
and  have  structural  defects;  these  characteristics  limit  bulk  mechanical 
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strength  and  modulus  rather  than  the  inherent  chain  properties. 

The  specific  properties  for  which  comparisons  will  be  made  are 
summarized  in  Table  1.  Each  property  provides  a  different  level  of 
comparison  with  the  SE  method.  The  equilibrium  geometry  is  one  of  the 
four  properties  against  which  the  SE  method  was  compared  and  gives  an 
indication  of  the  baseline  error  involved.  Because  polymer  modulus  is  the 
net  result  of  many  deformations  in  the  repeat  unit  (bond  extension,  bond 
angle  bending,  torsion  angle  twisting),  modulus  alone  is  not  a  sensitive  test 
of  the  SE  method.  In  order  to  determine  the  effect  of  the  SE  method  on  the 
calculated  modulus,  one  must  compare  SE  predictions  for  specific 
components  within  the  repeat  unit  with  experiment.  In  order  to  obtain 
information  on  specific  bonds  within  the  molecule,  one  needs  to  examine 
spectroscopic  data,  the  energy  shift  of  Raman  and  IR  active  lines  as  a 
function  of  strain.  Spectroscopy  provides  information  on  second  and  third 
derivatives  of  bond  potentials  generated  by  the  SE  paramaterization,  which 
describe  the  bond  strength.  In  this  study,  comparisons  will  concentrate  on 
backbone  vibrational  modes  of  PE  and  PDA  which  show  the  largest  strain 
sensitivity.  The  strain-dependent  intra-molecular  potentials  calculated  by 
the  SE  method  will  be  directly  compared  with  empirical  Morse  potentials 
that  have  been  shown  to  agree  with  the  spectroscopic  data  on  the  polymers. 
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TABLE  1 


SUMMARY  OF  COMPARISONS 


PROPERTY  EE 

Equilibrium  Geometry 

Crystallographic  Results  * 

Ab-Initio  Calculation  * 

Substituent  Effect 
Modulus 

Bulk  Measurements  * 

X-Ray  Simulation  * 

Substituent  Effect a 
Olig  ner  Comparison 
Spectroscopy 

Strain  Dependent  Shift  * 

Molecular  Potentials  * 

Morse  C-C  Comparison  * 

Simple  Molecule  C-C  * 

Empirical  Force  Constants  * 


PDA 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


a  Substituents  in  PDA  refer  to  the  fact  that  while  PDA  has  a  simple  carbon 
backbone,  it  is  synthesized  with  various  large  chemical  groups  bonded  to 
the  backbone.  Calculations  were  initially  performed  with  hydrogen  atoms 
bonded  to  the  backbone  to  make  the  calculations  simpler. 
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In  this  thesis  I  will  first  outline  the  Hartree-Fock  method  and  the 
additional  approximations  made  in  the  semi-empirical  approaches.  The  ab- 
initio  and  semi-empirical  methods  have  been  exhaustively  documented  in 
the  scientific  literature,  and  the  reader  will  be  directed  to  appropriate 
references  for  complete  treatments.  An  outline  of  the  approaches  is 
illuminating  because  it  lays  the  foundation  for  the  ultimate  applicability  of 
the  semi-empirical  method  to  the  question  at  hand.  The  accuracy  of  the  SE 
formulation  used  in  this  work  to  its  calibration  set  of  compounds  will  be 
discussed,  as  well  as  the  procedure  used  to  obtain  the  data  from  which  the 
polymer  mechanical  properties  is  obtained.  Next,  a  brief  description  of  the 
pertinent  mechanical  properties  and  the  various  techniques  by  which  they 
are  measured  will  follow,  as  well  as  the  calculation  of  these  properties  from 
the  SE  results.  The  results  of  the  calculations  on  PE  will  be  presented,  as 
well  as  comparisons  with  experimental  results.  The  results  of  the 
calculations  on  PDA  will  then  be  presented.  Following  this,  an  analysis  of 
the  results  and  a  discussion  of  the  implications  of  the  results  to  the  central 
question  of  this  thesis,  the  appropriateness  of  the  SE  technique  to  the 
calculation  of  polymer  mechanical  properties  away  from  equilibrium,  will 
be  presented.  These  comparisons  provide  evidence  that  SE  calculations 
systematically  overpredict  the  modulus  by  20%  to  25%. 
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II  APPROXIMATE  MOLECULAR  STRUCTURE  CALCULATIONS 
Ab  Initio  Method 

The  Schrodinger  wave  equation  for  atomic  hydrogen,  the  archetypal  one- 
electron  system,  describes  the  motion  of  the  electron  in  the  Coulomb 
potential  of  the  nucleus.  The  wave  equation  is  an  eigenvalue  problem  of  the 
form 

H  I  ¥  >  =  E  I  V  >  (1) 


where  H  is  the  Hamiltonian  operator.  For  atomic  hydrogen,  H  is  simply 


(2)  (Atomic  Units) 


where  2  the  kinetic  energy  operator  and  ^  is  the  spherically 
symmetric  Coulomb  potential  of  the  nucleus.  This  equation  is  separable  in 
spherical  polar  coordinates,  and  the  eigenfunctions  are  products  of 

spherical  harmonics  in  the  angular  variables  ©  and  0  and  associated 


Laguerre  polynomials  in  r.  The  wavefunctions  take  the  form 

H'n  1  m  (r,  0,  <}))  =  Rnl(r)  MW©)  eim0  (3). 


The  indices  arise  from  the  allowed  values  of  the  separation  constants,  and 
the  different  solutions  describe  allowed  discrete  states  for  the  hydrogen 
electron.3  The  wavefunction  corresponding  to  an  allowed  combination  of 
the  indices  n,l,m  is  referred  to  as  an  orbital  which  the  electron  is  said  to 
occupy.  One  follows  the  same  procedure  to  solve  the  wave  equation  for  a 
system  of  interacting  electrons  and  nuclei,  except  that  separable  solutions 
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cannot  be  obtained.  The  Hamiltonian  for  a  system  of  n  electrons  and  M 
nuclei  is 


ht=  M-  ivi 


n  M 


n  n 


**  rr  ..  ..  MM 

-  XX#^  +  XX^  + 


i=l  Z  A=1  Z  i=l  A=1  "lA  i=l  j>i 


rij 


A=1  B>A 


(4) 


n 


m 


IV 


v 


I  Kinetic  energy  of  the  electrons 

I I  Kinetic  energy  of  the  nuclei 

III  Coulomb  potential  between  electrons  and  nuclei 

I V  Electron-electron  repulsion 

V  Nuclei-nuclei  repulsion. 

The  wave  equation  for  the  system  is 

H-rotal  TO, 2  =  E-rotal  TO, 2  ...n;l,2,...M)  (5) 

where  'F  is  the  wavefunction  for  all  particles  in  the  molecule  and  E  is  the 
total  energy  of  the  system.  Since  each  particle  has  three  degrees  of 
freedom,  this  is  a  nonseparable  partial  differential  equation  in  3n  +  3M 
variables.4 

One  now  applies  several  approximations  to  simplify  this  problem.  The 
first  is  the  Born-Oppenheimer  approximation,  which  asserts  that  because 
of  the  large  mass  difference  of  the  electrons  and  the  nuclei,  the  nuclei  can 
be  considered  to  be  fixed  on  the  timescale  of  electronic  motion.  Thus,  one 
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can  initially  neglect  the  kinetic  energy  of  the  nuclei  and  consider  the 
nuclear  repulsion  a  constant  (which  does  not  affect  the  solution  of  the 
electronic  eigenvalue  problem).  One  first  solves  the  wave  equation  for  the 
electronic  Hamiltonian,  which  gives  the  electronic  wavefunctions  that 
depend  explicitly  on  the  electron  coordinates  but  only  parametrically  on  the 
nuclear  coordinates.  After  the  electronic  problem  has  been  solved,  the 
motion  of  the  nuclei  is  calculated  using  the  average  value  of  the  Coulomb 
potential  generated  by  the  electrons.  Solutions  to  the  nuclear  wave  equation 
describe  the  vibration,  rotation,  and  translation  of  the  molecule. 

The  next  step  in  the  approximate  solution  of  the  multi-electron  wave 
equation  is  to  seek  a  form  of  the  multi-electron  wavefunction  from  a 
combination  of  functions  which  depend  on  the  coordinates  of  only  one 
electron.  The  simplest  way  to  do  this  is  to  write  the  total  wave  function  as  a 
product  of  one  electron  wave  functions 

<m,2,3...n)  =  (!>i(l)  te(2)taP) ....  <j>N(n)  (6) 

The  probability  density  function,  T2,  which  describes  the  average  location  of 

the  electrons  is  then  a  product  of  the  one  electron  probabilities  <J>j 2.  This 

description  is  accurate  only  when  the  events  associated  with  each  of  the 
probabilities  occur  independently  of  each  other.  Thus,  the  physical  model 
used  in  this  approximation  is  that  the  electron  motions  are  independent  of 
each  other,  while  in  reality  correlation  of  electronic  motion  is  important.  A 
multi-electron  wavefunction  of  this  form  also  violates  the  anti-symmetry 
requirement  for  the  wavefunction  which  arises  from  the  Pauli  exclusion 
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prinriple.5  A  discussion  of  electron  spin  and  symmetry  requirements  of  the 
wave  functions  is  omitted  here;  these  considerations  require  more  complex 
combinations  of  one  electron  wave  functions  to  be  used  to  represent  the 
multi-electron  wavefunction.6 

Neglecting  electron-electron  repulsion  (in  addition  to  the  nuclear 
terms),  the  multi-electron  Hamiltonian  could  be  written  as  a  simple  sum  of 
one  electron  Hamiltonians.  Solutions  to  the  wave  equation  for  each  electron 
could  be  obtained  analytically  by  separation  of  variables,  and  the  exact 
solution  to  the  multi-electron  wavefunction  would  be  the  product  of  the  one- 
electron  wavefunctions.  Because  the  inter-electron  (pairwise)  repulsion 
depends  on  the  instantaneous  coordinates  of  two  electrons,  the  exact 
solution  cannot  be  expressed  in  a  separable  form.  However,  by  considering 
the  interaction  of  an  electron  i  with  all  the  other  n-1  electrons  as  an 
interaction  with  an  average,  effective  potential,  Fock7  was  able  to  derive 
approximate  one-electron  Hamiltonians  based  on  the  earlier  work  of 
Hartree8.  The  Hartree-Fock  equations  have  the  form: 

F  y=  [Hcore  +  (Ij  (2  Jj  -  Kj)]  'Fij  =  Xj  Eij  1=1, 2, 3,.. .n  (7) 

Here  F,  the  Fock  operator,  may  be  considered  an  effective  one-electron 
Hamiltonian  for  the  electron  in  the  molecular  environment,  and  its  various 
terms  have  simple  physical  interpretations.  Hcore  is  the  one-electron 
Hamiltonian  for  an  electron  moving  in  the  field  of  bare  nuclei.  The  term 
2Jj,  with  j  not  equal  i,  is  the  averaged  electrostatic  potential  of  the  two 

electrons  (of  opposite  spin)  in  the  orbital  Tj  .  The  exchange  potential  Kj 
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arises  from  the  indistinguishability  of  fermions. 

In  order  to  calculate  the  wavefunctions  using  the  Hartree-Fock 
equations,  one  needs  to  calculate  the  average  potential  from  n-1  electrons  in 
the  molecule,  which  is  determined  by  their  locations,  which  is  described  by 
their  wavefunctions,  i.e.,  you  need  to  know  the  answer  to  calculate  the 
answer.  The  general  process  for  solving  the  Hartree-Fock  equations  is  an 
iterative  process.  First,  one  assumes  a  set  of  trial  solutions,  Tj'.'lV,...  Tn' 

which  allows  computation  of  the  coulomb  and  exchange  operators  and 
thus  a  first  approximation  to  the  Fock  operator.  The  eigenfunctions  of  this 
operator  are  used  as  a  second  approximation  to  the  wavefunctions,  and  the 
procedure  is  repeated  until  the  Fock  operator  no  longer  changes  (to  a  given 
tolerance).  These  wavefunctions  are  said  to  be  self  consistent  with  the 
potential  field  they  generate,  so  this  procedure  is  known  as  the  self- 
consistent  field  (SCF)  method. 

The  choice  of  the  trial  wavefunctions  is  important  to  the  solution.  The 
most  commonly  used  choice  is  to  approximate  the  Hartree-Fock 
wavefunctions  with  linear  combinations  of  atomic  orbitals  (LCAO).  Since 
the  atomic  orbitals  (Eqn  3)  are  orthonormal  functions,  this  is  a  finite 
generalized  Fourier  series  expansion  of  the  molecular  orbitals.  This 
selection  has  the  further  advantage  in  aiding  the  interpretation  of  the 
results  since  the  molecular  properties  can  now  be  related  to  those  of  the 
constituent  atoms.  The  atomic  orbitals  used  to  construct  the  molecular 
orbitals  are  called  a  basis  set.  Better  representations  of  the  molecular 
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orbitals  can  be  obtained  by  including  more  terms  in  the  expansion,  a  larger 
basis  set,  though  this  increases  the  computation  time.  In  practice,  the 
atomic  orbitals  themselves  are  approximated  by  a  series  expansion  because 
the  actual  atomic  wavefunctions  are  computationally  difficult  to  evaluate. 
A  common  choice  is  to  use  Gaussian  functions;  strictly  speaking  this  is  not 
a  proper  choice  because  Gaussians  are  not  orthogonal  functions;  however, 
the  basis  set  can  be  orthogonalized  with  an  additional  operator. 

When  you  start  adding  up  all  the  terms  in  the  molecular  orbitals  it 
becomes  apparent  why  ab-initio  Hartree-Fock  calculations  become  complex 
very  quickly.  For  example,  consider  the  diatomic  molecule  HF.  One  needs 
six  atomic  orbitals  to  represent  each  molecular  orbital  .  If  one  uses  6 
Gaussian  functions  to  represent  each  atomic  orbital,  one  would  have  36  X 
36  X  6  complicated  integrations  to  perform  on  each  iteration  of  the  SCF 
procedure.  In  practice,  many  of  these  integrals  do  not  need  to  be 
recalculated  on  each  iteration, and  can  be  stored  on  disk;  however,  a  high 
level  calculation  on  a  fairly  simple  system  can  easily  consume  gigabytes  of 
disk  space.  There  are  many  other  considerations  that  make  the  ab-initio 
implementation  of  the  Hartree-Fock  method  computationally  expensive, 
both  in  terms  of  CPU  time  and  disk  space,  which  are  omitted  from  this 
discussion  for  the  sake  of  brevity.  The  net  result  is  that  while  ab-initio 
calculations  give  very  good  results  (when  compared  to  experiment),  the 
computational  cost  is  prohibitive  for  all  but  fairly  small  molecules. 

It  should  be  noted  that  because  the  Hartree-Fock  method  explicitly 
neglects  electron  correlation,  the  result  of  a  Hartree-Fock  calculation,  even 
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with  a  basis  set  of  infinite  size,  will  always  be  higher  than  the  true 
minimum  of  the  system.  Inclusion  of  electron  correlation  tends  to  decrease 
the  calculated  bond  stiffness,  which  implies  that  even  ab-initio  Hartree- 
Fock  calculations  of  mechanical  properties  may  be  systematically  high.9.10 
A  complete  discussion  of  the  theory  of  ab-initio  Hartree-Fock  method  is 
found  in  Szabo  and  Ostlund.6 

Semi-Empirical  Method 

In  order  to  enable  the  calculation  of  larger  molecules,  additional 
approximations  are  required  than  those  used  in  the  ab-initio  method.  The 
most  difficult  and  time  consuming  procedure  of  the  ab-initio  method  is  the 
computation  of  the  large  number  of  two  electron  (repulsion  or  exchange) 
integrals.  One  can  argue  from  physical  considerations  that  many  of  these 
integrals  should  be  very  small,  a  fact  which  can  be  verified  by  doing  the 
calculation.  One  important  class  of  electron-electron  interactions  that  is 
small  is  the  overlap  between  different  atomic  orbitals  describing  the 
molecular  orbital  of  an  electron,  e.g.,  the  overlap  between  a  IS  and  2P 
orbitals  from  the  same  atom  is  essentially  zero.  The  overlap  between  two 
different  atomic  orbitals  describing  the  molecular  orbital  of  an  electron  is 
referred  to  as  differential  overlap  (DO).  If  the  atomic  orbitals  originate  from 
the  same  atom,  it  is  called  monatomic  differential  overlap:  if  from  different 
atoms  it  is  called  diatomic  differential  overlap.  By  assuming  differential 
overlap  to  be  negligible,  a  large  number  of  integrals  involved  in  the  Fock 
operator  can  be  set  to  zero  (Neglect  of  Differential  Overlap).  This 


assumption  formed  the  basis  for  the  first  semi-empirical  method,  the 
complete  neglect  of  differential  overlap  (CNDO)  method  introduced  by  Pople, 
Santry,  and  Segaln.  In  addition  to  NDO,  only  valence  electrons  are  treated 
explicitly,  while  inner  shell  electrons  are  treated  as  part  of  the  nuclear 
core,  modifying  the  nuclear  potential  in  the  core  Hamiltonian  (Hc°re). 
Additional  based  semi-empirical  methods  include  intermediate  neglect  of 
DO  (INDO),  modified  intermediate  neglect  of  DO  (MINDO),  modified 
neglect  of  DO  (MNDO),  etc.  11 

Another  essential  simplification  in  the  SE  method  is  the  treatment  of 
interaction  integrals  in  the  Hamiltonian  as  element  specific  parameters  . 
The  number  of  parameters  required  to  describe  an  element  corresponds  to 
the  number  of  interaction  terms  in  the  approximate  Hamiltonian  used. 
Examples  of  parameters  would  be  the  one-electron  energy  of  an  atomic 
orbital  of  an  ion  (bare  nucleus  +  core  electrons)  resulting  from  the  removal 
of  all  valence  electrons,  an  atomic  Slater  orbital  exponent,  or  the  exponent 
in  a  gaussian  describing  core-core  repulsion.  The  name  semi-empirical 
arises  because  some  of  these  parameters  are  obtained  from  experimental 
data. 

Within  the  SE  approach,  the  self-consistent  field  algorithm  is  used  to 
calculate  the  molecular  orbitals  while  the  calculation  of  a  large  number  of 
complicated  integrals  in  the  ab-initio  method  is  replaced  by  a  relatively 
small  lookup  table.  Thus  a  semi-empirical  calculation  is  simplified,  as  is 
shown  in  Fig  1,  by  a  comparison  of  the  ab-initio  (la)  and  SE  (lb)  approaches 
which  is  reproduced  from  Clarkiz.  The  combination  of  greatly  reduced 
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(a)  (b) 

Figure  1.  Typical  flow  charts  for  (a)  an  ab-initio  optimization  and  (b)  a 
corresponding  semi-empirical  calculation. 
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computation  times  and  memory  requirements  means  that  much  larger 
systems  can  be  treated  by  SF  methods  than  by  ab-initio  techniques. 
Moreover,  the  accuracy  of  the  SE  methods  in  calculating  heats  of  formation 
compares  favorably  on  a  wide  range  of  systems  when  compared  to  moderate 
level  ab-initio  calculations.  Note,  however,  that  this  is  not  true  for  all 
molecular  properties.  13 

The  ultimate  success  of  a  semi-empirical  method  hinges  on  the  validity 
of  the  approximate  Hamiltonian  used  and  on  the  values  of  the  parameters 
used.  Ideally,  one  would  want  the  parameters  that  describe  each  element 
to  be  completely  independent  of  the  molecular  environment,  i.e.  the 
parameters  which  describe  a  carbon  atom  in  CH4  should  give  results  of  the 

same  accuracy  when  used  in  CF4.  In  reality,  the  specific  molecular 

environment  does  affect  the  atomic  parameters.  In  order  to  make  the 
application  of  the  method  generic  the  parameters  are  selected  to  minimize 
the  least  squares  error  for  selected  molecular  properties  over  a  large  set  of 
molecules. 

It  should  be  noted  that  the  use  of  experimental  data  to  define  some  of  the 
parameters  actually  gives  the  semi-empirical  method  some  advantage  over 
the  ab-initio  method.  The  measured  quantities  inherently  contain  all  the 
physics  involved  and  so  are  not  limited  by  the  approximations  used  to  define 
a  tractable  Hamiltonian.  A  complete  discussion  of  SE  methods  is  presented 
by  Pople  and  Beveridge.  4 

The  Austin  Model  1  (AMI)  semi-empirical  Hamiltoniar  (and 
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parameters)  were  used  in  the  calculations  in  this  thesis.  AMI  is  an 
extension  of  the  earlier  MNDO/3  Hamiltonian  which  is  known  to 
overestimate  nuclear  core-core  repulsion.  Thus  the  parameters  for  the 
AMI  include  those  from  the  approximate  Hamiltonian  as  well  as 
Gaussian  exponents  describing  the  modified  core-core  repulsion.  Five  of  the 
parameters  were  assigned  values  from  atomic  spectroscopy.  The  other 
parameters  were  initially  estimated  and  adjusted  to  best  reproduce  the 
following  experimental  equilibrium  (297  K)  gas  phase  properties: 

1.  Heat  of  Formation 

2.  Dipole  Moments 

3.  Ionization  Potentials 

4.  Molecular  Geometries. 

Most  of  these  properties  used  in  the  parameterization  are  obtained  from 
experiment,  though  for  some  molecules  on  which  precise  measurements 
are  impractical,  the  results  of  high  level  ab-initio  calculations  are 
considered  reliable  enough  to  be  used.  Table  2,  which  reproduces  Table  XII 
from  Stewart  14  shows  the  average  error  in  the  calculation  of  molecular 
geometries  associated  with  the  AMI  method.  A  detailed  statistical 
analysis  of  the  errors  associated  with  AMI  with  respect  to  these  four 
properties  is  presented  in  ref  11.  Of  these  benchmark  properties  only 
comparisons  with  equilibrium  molecular  geometry  will  be  performed  in 
this  study. 
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TABLE  2 


AVERAGE  ERRORS  IN  MOLECULAR  GEOMETRIES 


Geometric 

Parameter 

Number  of 
Molecules  for 
Comparison 

MNDO 

AMI 

PM3* 

Bond  Length  [A] 

372 

0.054 

0.050 

0.036 

Bond  Angle  [Deg] 

158 

4.342 

3.281 

3.932 

Torsion  Angle  [Deg] 

16 

21.619 

12.494 

14.875 

*The  PM3  Hamiltonian  is  a  modification  of  MNDO  which  gives  better 
results  for  hypervalent  compounds. 

Semi-Empirical  Calculations  of  the  Electronic  Structure  of  Polymers 

While  large  molecules  can  be  handled  by  the  semi-empirical  method,  an 
ideal  polymer  is  an  infinite  chain  molecule;  clearly  an  approximation  must 
be  made  to  calculate  polymer  properties,  even  within  the  SE  approximation. 
One  technique  is  to  consider  a  series  of  small  molecules  which  are  made 
up  of  increasing  numbers  of  the  polymer  repeat  units,  namely  oligomers. 
The  oligomer  properties  will  not  accurately  represent  the  polymer 
properties,  but  one  can  extrapolate  the  results  as  a  function  of  repeat  units 
to  the  infinite  repeat  unit  result.  Drawbacks  to  this  method  are  that  it  is 
time  consuming  and  end  effects  will  always  be  present  in  the  oligomers 
which  are  not  in  found  in  polymers.  Another  approach  is  the  cluster 
method  developed  by  Perkins  and  Stewart.15  In  the  cluster  method,  one 
performs  calculation  on  a  cluster,  which  is  a  number  of  polymer  repeat 
units.  The  atom  at  one  end  of  the  cluster  interacts  with  the  atom  at  the 


other  end,  as  though  it  were  connected  in  a  ring.  Ref  15  provides  a  good 
discussion  of  how  this  approach  for  a  cluster  of  carbon  atoms  reduces  to 
Huckel  theory  for  benzene.  Using  the  cluster,  exact  results  for  an  infinite 
molecule  can  be  reproduced  with  as  few  as  one  repeat  unit,  provided  the 
repeat  unit  is  long  enough  to  make  an  acceptable  cluster.  The  number  of 
repeat  units  required  in  the  cluster  is  driven  by  the  electron  de-localization 
length  in  the  polymer;  the  length  of  the  cluster  must  be  large  enough  so 
that  any  one  atom  does  not  interact  with  itself  through  the  cyclic  boundary 
conditions.  Stewart  suggests  that  a  cluster  length  of  10  A  is  sufficient  for 
most  polymers,  except  those  containing  conjugated  Pi  bonds,  where  a 
cluster  length  of  20  A  is  recommended.  This  was  verified  in  this  study  for 
PDA,  where  the  heats  of  formation  and  repeat  unit  lengths  were  monitored 
as  a  function  of  cluster  size.  PE  has  been  previously  calculated  using  the 
cluster  method  with  the  MNDO/3  Hamiltonian,  so  that  the  cluster  length 
dependence  is  known.  1 
Details  of  the  Calculations 

The  program  MOPAC  version  5.0  was  used  with  the  AMI  Hamiltonian 
was  used  in  this  study  with  the  unrestricted  Hartee-Fock  (UHF)  method.  16 
The  cluster  method  as  described  in  ref  15  was  used  to  simulate  ideal 
polymer  chains.  Strain  dependent  geometry  optimizations  were  performed 
by  using  the  translation  vector,  which  describes  the  length  of  the  cluster 
and  defines  the  polymer  connectivity,  as  a  reaction  coordinate  for  a  series  of 
strain  values.  All  MOPAC  calculations  were  performed  on  a  Silicon 
Graphics  workstation  in  the  Polymer  Branch  of  the  Materials  Laboratory  at 
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WPAFB.  Ab-initio  oligomer  calculations  were  performed  using  the 

programs  GAUSSIAN88  and  GAUSSIANSO1?  on  Elxsi  and  Cray  supercomputers. 
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Ill 


RELATIONSHIP  OF  THEORY  TO  EXPERIMENT 


Calculations  of  Mechanical  Properties 

The  basic  method  used  to  determine  the  mechanical  properties  of  a 
polymer  is  to  measure  the  amount  of  force  required  to  change  the  length  of 
a  sample.  From  this  measurement  one  obtains  a  stress  (force  per  unit 
area)  strain  (fractional  change  in  length  away  from  equilibrium)  curve.  A 
generic  stress-strain  curve  is  shown  in  Fig  2.  Several  mechanical 
properties,  defined  by  the  behavior  of  the  stress-strain  curve,  are  labelled. 
The  modulus  is  the  initial  slope  of  the  stress-strain  curve.  The  yield 
strength  is  the  amount  of  force  at  the  point  where  the  slope  of  the  stress- 
strain  goes  to  zero.i.e.,  the  resistance  to  deformation  vanishes.  The 
ultimate  strength  is  the  force  at  the  point  where  the  material  physically 
breaks.  If  one  had  a  perfectly  ordered,  perfectly  aligned,  defect-free 
crystalline  polymer,  the  stress-strain  curve  measured  on  the  bulk  material 
would  be  identical  to  the  properties  of  a  single  polymer  chain. 

Straining  a  moleculat  system,  i.e.,  increasing  or  decreasing  the 
separation  between  the  atoms,  requires  energy,  so  that  the  potential  energy 
of  the  system  increases  when  it  is  constrained  to  a  non-equilibrium  length. 
To  a  first  approximation,  the  interaction  between  two  atoms  in  a  molecule 
can  be  described  as  a  harmonic  oscillator  (ideal  spring).  For  an  ideal 
harmonic  oscillator  (H.O.),  this  energy  is  proportional  to  the  square  of  the 
deviation  from  the  equilibrium  separation.  The  force  required  to  stretch  (or 
squeeze)  a  spring  is  proportional  to  the  second  derivative  of  the  potential 
energy  curve  with  respect  to  position.  The  second  derivative  of  a  quadratic 
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Figure  2.  Generalized  tensile  stress-strain  curve  for  a  polymer 
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potential  is  a  constant,  so  the  stiffness  of  a  spring  is  described  by  a  force 
constant.  The  resistance  of  a  material  to  deformation  is  described  in  terms 
of  its  modulus,  the  force  per  unit  area  required  to  distort  it  from  its 
equilibrium  length  (along  a  specific  axis),  which  is  equal  to  the  initial  slope 
of  the  materials'  stress-strain  curve.  The  modulus  of  a  spring  then  is  the 
second  derivative  of  its  energy  potential  (force  constant)  times  its 
equilibrium  length  divided  by  its  cross  sectional  area. 

The  energy  potential  curve  of  a  strained  molecular  system  can  be 
generated  through  semi-empirical  calculations.  First,  one  determines  the 
equilibrium  geometry  of  the  molecule  by  optimizing  all  its  variables  (bond 
lengths,  bond  angles,  torsion  angles).  Associated  with  this  optimum 
equilibrium  geometry  is  a  heat  of  formation.  One  then  changes  a  variable 
which  describes  the  overall  molecular  length,  e.g.  increase  the  length  by 
2%,  and  then  calculates  the  optimum  geometry  for  the  molecule  with  this 
variable  fixed.  Clearly,  the  heat  of  formation  for  this  strained  geometry  will 
be  higher  than  equilibrium.  By  calculating  the  heats  of  formation  for  a 
number  of  such  strain  values  one  can  map  out  the  potential  energy  curve 
for  the  molecule.  The  second  derivative  of  this  curve  is  the  spring  constant 
of  the  molecule,  which  can  be  used  to  calculate  its  modulus  (on  a  molecular 
level,  the  crystal  unit  cell  area  is  used  for  the  cross  sectional  area).  For  a 
polymer,  one  fixes  the  length  of  the  cluster,  optimizes  the  geometry  of  the 
repeat  unit  and  obtains  the  minimum  energy  for  a  particular  strain.  For  a 
polymer,  this  energy  is  an  amalgam  of  all  the  various  deformation  modes 
(bond  extensions,  bond  angle  changes,  and  torsion  angle  twists). 


Comparing  the  calculated  cluster  modulus  with  experiment  is  not  a 
sufficient  test  of  the  validity  of  the  calculation  because  compensating  errors 
in  the  calculated  deformation  modes  could  combine  to  produce  a  value  for 
the  modulus  which  seems  "reasonable",  i.e.,  is  less  than  several  orders  of 
magnitude  larger  than  experiment. 

There  is  a  potential  problem  inherent  with  this  technique:  the 
parameters  used  in  AMI  were  optimized  to  reproduce  the  equilibrium  heat 
of  formation.  There  is  no  guarantee  that  the  heats  of  formation  are 
meaningful  away  from  equilibrium.  One  may  have  an  accurate  value  of 
the  potential  at  equilibrium,  but  the  optimization  procedure  includes 
information  about  the  quality  of  the  energy  derivatives.  Because  the 
mechanical  properties  of  a  material  depend  on  the  derivatives  of  the 
potential,  in  order  to  assess  the  accuracy  of  the  calculation  of  a  material's 
mechanical  properties,  one  must  investigate  the  accuracy  of  the  potentials 
(and  their  derivatives)  used  to  describe  the  molecular  bonding  in  the 
material. 

A  more  precise  description  of  a  molecular  bond  is  that  only  for  small 
deviations  about  the  equilibrium  length,  the  potential  is  harmonic  with 
displacement.  At  larger  displacements  the  true  anharmonicity  of  the 
potential  can  be  observed.  In  general,  the  potential  is  softer  than  harmonic 
in  tension  (positive  strain)  and  harder  in  compression  (negative  strain). 
Thus  the  spring  "constant"  for  a  real  bond  is  not  a  constant.  Common 
attempts  to  represent  the  anharmonicity  in  a  real  bond  potential  are  to 
include  cubic  or  higher  polynomial  terms.  A  very  good,  generic  fit  to  bond 


potentials  is  given  by  the  Morse  potential,  which  was  originally  derived  as 
an  empirical  fit  to  diatomic  spectra.  The  form  of  the  Morse  potential  is*8 

V(r)  =  D0  (1  -  e-a(r'ro))2  (8) 

where  Do  is  the  dissociation  energy  of  the  bond,  ro  is  the  equilibrium 

length,  and  a  describes  the  curvature  of  the  potential.  A  comparison  of  a 
Morse  and  harmonic  potentials  is  shown  in  Fig  3. 

The  vibrational  frequency  of  a  diatomic  molecule  is  given  by 


where  k  is  the  force  constant  and  p  is  the  reduced  mass  of  the  two  atoms. 

For  a  harmonic  potential,  the  value  of  k  does  not  change  so  the  frequency 
does  not  change  with  inter-atomic  separation.  For  a  system  described  by  an 
anharmonic  potential,  the  frequency  changes  with  strain  because  the 
second  derivative  of  the  pote*-  ial  is  no  longer  constant.  By  measuring  the 
vibrational  spectrum  of  a  molecule  one  obtains  the  value  of  the  second 
derivative  of  the  potential  at  equilibrium.  Measuring  the  shift  of  the 
vibrational  frequencies  with  strain  provides  information  about  the 
anharmonicity  of  the  potential.  Thus  by  comparing  measured  vibrational 
shifts  with  those  calculated  by  AMI  it  is  possible  to  determine  the  accuracy 
of  the  calculated  bond  stiffness. 

One  can  directly  calculate  the  strain  dependent  potential  of  a  single  bond 
in  a  molecule  by  constraining  just  that  particular  bond  length  and 
optimizing  the  rest  of  the  geometry.  This  proved  to  be  very  useful  for 
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Figure  3.  Comparison  of  Morse  and  harmonic  potentials.  Both 
have  the  same  second  derivative  at  equilibrium. 
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polyethylene  because  Wool19  provided  a  Morse  potential  for  the  C-C  bond 
stretch  which  was  used  to  calculate  the  modulus.  This  Morse  potential  is 
judged  to  be  a  good  representation  of  the  C-C  bond  because  it  matched  both 
the  absolute  frequencies  and  the  strain  dependent  frequency  shifts  observed 
in  polyethylene. 

It  is  known  that  the  calculation  of  vibrational  frequencies  from  the  SE 
(or  an  ab-initio)  method  produces  results  which  systematically  overestimate 
the  measured  frequencies.6  This  discrepancy  is  a  result  of  the  method  used 
to  calculate  the  frequency  rather  than  inherent  in  the  SE  method  itself.  The 
force  constant  is  calculated  by  taking  the  second  derivative  of  the  potential 
at  equilibrium.  The  vibrational  energy  levels  of  a  harmonic  oscillator  are 

E(V)  =  'Do  +  2^‘°(V  +  2)  (9) 

where  Do  is  the  minimum  of  the  potential,  v  is  a  non-negative  integer,  and 


to  = 


1  d2V(r) 


dr2  r=ro  (10) 

In  equation  (10  ),  |X  is  the  reduced  mass  of  the  molecule.  The  vibrational 
energy  levels  of  a  Morse  potential  are 


E(v)  =  -D0  +  -t-  cg{(v  + 1)  -  1<v  +  i)2} 
2k  2  c,  2 


(11) 


£  is  defined  such  that  the  anharmonic  (Cu  +  \)i)  term  can  never  exceed  the 


harmonic  ((o  +2  ))  term.20  It  is  seen  by  a  comparison  of  Eqns  (9)  and  (11) 
that  the  energies  of  the  harmonic  and  Morse  potentials  will  differ 
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significantly  only  if  u  is  large.  This  is  not  the  case  for  the  v=l  state  as 
calculated  by  the  SE  codes,  or  if  the  energy  separation  between  the 
minimum  of  the  potential  and  the  u=l  state  is  large,  which  occurs  if  p  is 

small.  This  approximation  in  the  frequency  calculation,  as  well  as  the 
neglect  of  electron  correlation,  produces  systematic  errors;  because 
comparisons  are  made  to  strain  dependent  frequency  differences,  the 
systematic  error  cancels  out  and  comparisons  with  experimental  strain 
dependent  frequency  shifts  is  meaningful. 

Matching  vibrational  frequencies  and  shifts  is  more  complicated  for 
polyatomic  molecules  than  diatomics  because  the  normal  modes  of 
vibration  involve  displacement  of  several  atoms  simultaneous1  y  (and 
because  the  number  of  vibrational  frequencies  scales  as  3n-6),  hence  it  is 
more  difficult  to  isolate  the  individual  bond  force  constants  involved. 
Although  an  ideal  polymer  has  an  infinite  number  of  atoms,  its  repetitive 
structure  results  in  many  degenerate  vibrational  modes,  and  so 
interpretation  of  the  spectra  is  tractable. 

Measurement  of  Polymer  Modulus 

As  stated  before,  the  measurement  of  the  bulk  mechanical  properties  for 
a  polymer  by  obtaining  its  stress-strain  curve  does  not  represent  the 
ultimate  mechanical  properties  for  that  polymer.  Fig  4  shows  the 
structure  of  a  typical  amorphous  polymer,  in  which  the  polymer  chains  are 
randomly  oriented  and  twisted  around  each  other.  When  one  strains  a  bulk 
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sample  of  an  amorphous  polymer,  one  is  addressing  primarily  the 
strengths  of  the  inter-chain  interactions.  Only  those  regions  whose  chain 
axis  is  aligru  i  with  the  axis  of  strain  will  resist  with  the  chain  strength. 

As  t  re  degree  of  order  is  increased  in  the  polymer,  by  better  aligning  the 
chain  along  one  axis,  the  mechanical  modulus  of  the  polymer  will 
inc  ease.  If  the  individual  chains  were  perfectly  aligned,  the  mechanical 
modulus  would  equal  the  ultimate  axial  modulus  (the  transverse  modulus, 
however,  would  be  minimized).  In  practi  ,  "perfect"  alignment  can  only 
be  obtained  in  a  limited  number  of  polymers,  of  which  PDA  is  one  which 
can  be  polymerized  in  single  crystal  form.  Mechanical  properties  in  PDA 
are  not  constrained  by  alignment  but  by  crystal  defects,  chain  ends  and 
other  chemical  imperfections.^1, 22 

The  calculations  of  polymer  moduli  are  performed  on  an  ideal,  infinite 
single  polymer  chaiu,  which  simulates  perfect  polymer  morphology  and 
neglects  the  role  of  interchain  forces.  While  this  scenario  is  unrealizable  in 
a  practical  sense,  it  is  useful  because  it  defines  the  upper  bound  for  the 
mechanical  properties  of  a  given  polymer.  Comparisons  of  this  upper  limit 
for  different  polymers  is  useful  because  if  the  degree  of  alignment  one  can 
obtain  in  processing  is  roughly  equal  for  all  polymers  (or  at  least  a  class  of 
chemically  similar  polymers),  the  polymer  with  the  stiffest  /strongest  chain 
will  make  the  stiffest/strongest  bulk  material.  An  analysis  of  the 
deformation  modes  can  also  give  insight  into  what  chemical  structures 
contribute  most  to  mechanical  strength  and  stiffness. 

Since  mechanical  measurement  docs  not  represent  the  ultimate  limits 
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of  a  polymer,  other  techniques  are  employed  to  make  molecular  level 
measurements.  X-ray  diffraction  is  used  to  monitor  the  location  of  the  29 

peak  which  corresponds  to  the  chain  length;  in  PE  this  is  the  (002)  peak.  (In 
this  study,  the  polymer  axis  will  be  defined  as  the  c  crystal  axis.)  A 
calibrated  force  is  applied  to  the  polymer  to  elongate  it,  and  the  strain  of  the 
c-axis  is  monitored  from  the  appropriate  diffraction  peak.  Under  the 
assumption  that  stress  is  uniformly  distributed  throughout  the  sample 
(which  is  not  always  a  very  good  assumption)  one  obtains  a  molecular 
stress-strain  curve,  whose  slope  is  the  modulus.  Since  this  technique 
monitors  molecular  strain  rather  than  bulk  strain,  the  modulus  measured 
better  represents  the  ultimate  modulus  of  the  polymer.23  Bulk  PE  does  not 
have  a  pure  crystalline  structure,  so  X-ray  modulus  measurements  on  PE 
are  very  sensitive  to  the  degree  of  crystallinity.  The  X-ray  modulus  of  PDA 
should  equal  the  ultimate  theoretical  modulus  because  it  can  be 
polymerized  in  single  crystal,  though  defects  may  play  a  role  in  the 
measurement.  The  X-ray  modulus  of  PDA  has  not  been  reported  to  date. 

Another  method  of  measuring  the  ultimate  moduli  of  polymers  is  to 
derive  the  force  constants  associated  with  the  bonds  from  the  polymer's 
spectra  and  calculate  the  effective  axial  deformation  force  constant  from 
them.  This  technique  is  difficult  to  apply  in  practice  because  it  is  difficult  to 
perform  the  normal  vibrational  mode  analysis  on  the  polymer  and 
deconvolve  the  individual  force  constants  from  the  complex  spectra.  While 
this  method  has  no  predictive  capability  on  proposed  polymers,  in  theory  a 
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force  constant  analysis  should  give  the  correct  value  for  the  modulus 
because  it  is  based  on  measured  force  constants,  i.e..  there  are  no 
approximations  made  to  obtain  the  result.  In  practice,  for  all  but  simple 
polymers  the  mode  structure  is  so  complex  that  researchers  often  used 
measured  values  as  a  first  guess,  and  then  adjusted  the  force  constants  to 
fit  the  spectra.  The  set  of  force  constants  obtained  in  the  fitting  in  many 
cases  is  not  unique,  and  so  the  advantage  of  being  based  on  measurement  is 
lost. 

This  discussion  of  mechanical  properties  is  intended  to  be  a  generic 
description.  The  mechanical  behavior  of  polymers  is  much  more 
complicated,  and  differs  greatly  from  polymer  to  polymer  for  very  material 
specific  reasons.  A  good  introduction  to  the  mechanical  behavior  of 
polymers  is  presented  by  Young.20 


IV 


POLYETHYLENE  RESULTS  AND  DISCUSSION 


Molecular  Structure  and  Modulus  Calculations 

Polyethylene  (PE)  chains  are  comprised  of  single  carbon-carbon  bonds 
in  a  planar  zigzag  structure,  which  is  one  of  the  simplest  structures 
possible  for  any  polymer.  Comparisons  between  the  molecular  structure 
results  obtained  by  the  SE  and  ab-initio  methods  were  performed  on  the  PE 
oligomer  shown  in  Fig  5.  The  cluster  used  in  the  SE  polymer  calculations 
is  identical  to  the  oligomer  shown  except  that  there  are  only  two  hydrogens 
bonded  to  Cl  and  C6,  the  periodic  boundary  conditions  establishing  a  carbon- 
carbon  bond  at  the  appropriate  distance  and  angle  between  Cl  and  C6.  The 
SE  carbon-carbon  bond  distances  are  approximately  0.04  angstroms  shorter 
than  the  ab-initio  distances.  The  C-C  bond  distance  in  PE,  obtained  from 
X-ray  diffraction,  is  1.53  A.24  The  SE  bond  angles  are  approximately  2 
degrees  smaller  than  the  ab-initio  results.  The  SE  values  for  bond  lengths 
and  angles  are  accurate  to  within  the  published  uncertainly  of  the  AMI 
Hamiltonian  listed  in  Table  2. 

In  both  the  SE  and  ab-initio  oligomer  results,  bond  lengths  and  angles 
are  seen  to  depend  on  the  location  in  the  molecule,  with  the  bond  lengths  at 
the  end  of  the  oligomer  showing  the  largest  variation  (end  effects).  In 
contrast,  the  SE  results  with  periodic  conditions  are  essentially  equal 
throughout  the  cluster,  as  anticipated  for  an  ideal  polymer  chain. 

The  strain  dependent  heat  of  formation  potential  curve  for  PE  is  shown 
in  Fig  6  from  -10%  to  +10%  strain.  This  plot  is  offset  by  the  equilibrium 
heat  of  formation  of  the  cluster,  -38.9105  kCals/mole.  The  cluster's 
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Figure  5  .  Comparison  of  PE  oligomer  and  cluster  geometry  calculated 
by  AMI  and  ab-inition  methods. 


Figure  6.  PE  cluster  strain  dependent  heat  of  formation.  The 
discontinuity  in  the  curve  occurs  when  the  cluster  bends  out 

of  its  equilibrium  planar  geometry. 
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equilibrium  length  is  7.488  A.  The  discontinuity  in  the  curve  in 
compression  occurs  when  PE  deviates  from  a  planar  structure,  i.e.,  the  PE 
chain  bends.  This  deformation  mode  does  not  directly  correlate  to  a  bulk 
failure  mode  in  the  polymer,  however,  it  indicates  that  PE  does  not  possess 
a  large,  inherent  compressive  strength25.  The  PE  chain  has  not  failed  in 
tension  up  to  10%  strain;  again  this  does  not  correlate  to  a  bulk  polymer 
property.  Since  PE  has  a  planar  chain  structure  in  the  calculation,  the  only 
tensile  failure  mode  accessible  in  the  calculation  is  bond  breaking.  This 
calculation  predicts  that  a  PE  chain  will  not  fail  at  a  uniformly  distributed 
tensile  strain  of  10%. 

The  second  derivative  of  the  strain-dependent  heat  of  formation  curve  is 
the  restoring  force  of  the  cluster.  This  can  be  obtained  several  ways.  One 
can  perform  a  polynomial  fit  to  the  potential;  then  twice  the  second  order 
coefficient  is  the  force  constant  at  equilibrium  (k).  The  chain  modulus  can 
then  be  defined  by 

kxL 

A  do) 

where  L  is  the  cluster  length  and  A  the  cross  sectional  area.  In  this  study, 
the  cross  sectional  area,  obtained  from  the  crystal  structure,  was  taken  to 
remain  constant  (18.4  A2  24)  as  a  function  of  strain  (Poisson  ratio  of  0.0V 
Klei  and  Stewart  used  a  third  order  polynomial,  which  is  the  first  order 
anharmonic  perturbation  to  a  harmonic  potential,  to  fit  the  potential.  The 
second  derivative  of  a  third  order  polynomial  is  a  line,  thus  a  cubic  fit 
predicts  that  the  restoring  force  drops  off  linearly  with  strain.  One  can  take 


the  second  derivative  of  the  potential  numerically  and  attempt  to  obtain  the 
exact  strain  dependent  behavior  of  the  modulus  predicted  by  the 
calculation.  While  the  potential  (Fig  6)  appears  smooth,  random  variations 
in  the  heat  of  formation  values  requires  that  smoothing  be  used  in 
conjunction  with  the  numerical  differentiation.  Typically,  the  first 
derivative  of  the  potential,  the  force  curve,  is  obtained  directly.  The  force 
curve  is  smoothed  with  a  fifth  order  Gaussian  filter  before  the  second 
derivative  is  taken.  The  second  derivative  curve  is  also  smoothed.  A 
comparison  of  the  numerical  second  derivative  and  the  second  derivative 
from  a  cubic  polynomial  fit  is  shown  in  Fig  7.  The  units  of  the  second 
derivative  have  been  converted  to  SI,  newtons  per  meter  (Nt/m),  in  this 
figure.  As  is  seen,  both  procedures  agree  very  closely  at  equilibrium,  with 
a  value  of  100  ±  5  Nt/m. 

The  modulus  of  PE  predicted  by  this  AMI  calculation  is  400  ±  20  GPa. 
Reference  1,  which  used  the  MNDO  semi-empirical  Hamiltonian, 
predicted  a  value  of  360  GPa.  This  difference  is  consistent  with  the 
difference  in  the  heats  of  formation  obtained  with  these  different 
Hamiltonians.  An  earlier,  low  level  ab-initio  calculation  of  PE  chain 
modulus  reported  a  value  of  approximately  400  GPa.26  The  reported  value 
for  PE  chain  modulus,  calculated  and  measured  by  the  techniques 
discussed  previously,  varies  from  180  to  400  GPa.  27  The  commonly  accepted 
value  of  PE  chain  modulus  is  around  300  GPa.  While  the  AMI  value  of  400 
GPa  is  high,  it  is  not  outside  the  range  of  previously  reported  values. 
However,  a  priori  there  is  no  justification  to  claim  that  the  AMI  result  is 
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Figure  7  Comparison  of  second  derivatives  of  PE  heat  of 
formation  potential.  The  discontinuity  in  the  numerical  derivative 
results  from  molecular  failure. 
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more  or  less  accurate  than  any  of  the  others. 

In  order  to  investigate  the  accuracy  of  the  calculation,  an  obvious 
question  is  to  inspect  the  strain  dependent  minimum  energy  geometries 
and  how  the  c-axis  strain  of  the  cluster  is  distributed  among  the  various 
components.  The  length  of  the  cluster  can  increase  by  increasing  the  bond 
lengths  and  bond  angle;  one  would  expect  to  observe  the  length  of  the 
cluster  increase  by  some  combination  proportional  to  the  resistance  of  each 
deformation  mode.  (The  discussion  here  is  limited  to  tension  because  the 
experiments  on  PE  used  for  comparison  did  not  report  any  results  in 
compression.)  A  commonly  used  rule  of  thumb  puts  a  ratio  of  100:10:1  in 
terms  of  the  amount  of  force  required  to  change  a  bond  length,  bond  angle, 
and  torsion  angle  respectively.  21  One  would  then  expect  polyethylene  to 
deform  primarily  by  changing  the  bond  angles,  with  a  smaller  contribution 
coming  from  increasing  the  bond  lengths.  (Since  the  equilibrium  structure 
of  PE  is  planar,  a  change  in  the  torsion  angle  is  considered  a  failure  mode.) 

Figure  8  shows  the  relative  variation  of  the  bond  angles  and  bond 
lengths  in  PE  under  tension  as  a  function  of  cluster  c-axis  strain.  The  bond 
angle  variation  is  seen  to  be  larger  than  the  bond  length  variation  ,  but  not 
by  a  factor  of  10.  The  average  angular  contribution  to  c-axis  strain  is 
obtained  by  performing  a  linear  fit  to  this  data,  is  0.51  ±0.05,  while  the  bond 
length  contribution  is  0.46  ±  0.05.  This  shows  that  the  c-axis  deformation  is 
roughly  equally  distributed  between  these  two  modes,  not  a  10  to  1 
partitioning.  This  behavior  was  also  reported  by  Wool  17,  and  is  consistent 
with  observed  spectral  line  shifts. 
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\  Figure  8.  PE  distribution  of  chain  strain  between  bond  angle 

p  and  bond  length. 
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Figure  9.  PE  cluster  strain  potentials  with  bond  angle  and  bond  length 
fixed.  The  unconstrained  potential,  i.e.,  nothing  fixed,  is  shown  for 
comparison. 


42 


As  a  check  on  this  behavior,  the  modulus  for  PE  was  recalculated  for 
two  constrained  cases:  one  in  which  the  angles  were  held  fixed  so  the  c-axis 
strain  could  only  come  from  changes  in  the  bond  length,  and  one  in  which 
the  bond  lengths  were  held  fixed  so  that  c-axis  strain  could  only  come  from 
changes  in  the  bond  angles.  The  potential  curves  for  these  two  cases  are 
shown  in  Fig  9.  The  difference  between  these  two  potentials  is  seen  to  be 
small.  The  bond  extension  modulus  obtained  is  525  GPa,  and  the  bond 
angle  modulus  is  508  GPa.  The  difference  between  these  two  values  is 
within  the  accuracy  of  the  calculation.  This  result  is  consistent  with  the 
previous  one,  that  there  is  a  small  difference  between  bond  angle 
deformation  and  bond  length  extension.  However,  both  restricted 
deformations  are  stiffer  than  the  fully  free  deformation,  520  GPa  to  400  GPa. 
The  substantial  difference  in  the  constrained  and  unconstrained  moduli 
can  be  explained  by  the  fact  that  bond  lengths  and  bond  angles  are  not 
independent  variables  in  a  quantum  mechanics  calculation,  but  are 
manifestations  of  the  minimum  energy  electronic  structure  of  a  molecule. 

Simulations  of  the  X-ray  diffraction  patterns  were  also  performed  as  a 
check  on  the  strain  dependent  PE  geometries.  The  crystal  structure  of  PE 
is  shown  in  Fig  10.  Crystal  PE  is  orthorhombic,  Pnam  symmetry,  with  unit 

cell  dimensions  of  a=7.417  A,  6=4.954  A,  and  c=2.534  A24.  Because  there 
are  two  chains  per  unit  cell,  the  cross  sectional  area  is  (a  x  6)  /  2  .  The 

location  of  the  (006)  20  reflection  peak  of  the  C6H12  cluster  describes  plane 
spacing  only  along  the  chain  axis;  becsause  an  infinite  number  of 
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PE  crystal  structure  .  Here  c  ic  used  as  the  polymer  axis. 
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Figure  11.  Location  of  006  X-ray  diffraction  peak  as  a  function  of  c  -axis  strain 
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combinations  of  bond  angles  and  bond  lengths  can  produce  the  same 
spacing  along  the  chain  axis,  the  (006)  peak  location  alone  cannot  be  used  to 
infer  information  about  PE  deformation.  The  location  of  the  (006)  reflection 
as  a  function  of  strain  is  shown  in  Fig  11.  The  modulus  can  be  calculated 
from  the  relation 

e  =  Astress. 

M 

d  (12) 

where 


d  _  X  (X-ray) 

2  sm(  (13) 

The  X-ray  modulus  calculated  for  PE  is  400  GPa,  which  is  consistent  with 
the  modulus  calculated  from  the  second  derivative  of  the  AMI  potential. 
Since  the  first  derivative  of  this  potential  is  used  to  define  the  stress  curve, 
the  X-ray  simulation  is  really  just  a  consistency  check,  not  an  independent 
calculation  of  the  modulus. 

Comparisons  with  Spectroscopy 

The  calculation  of  the  vibrational  spectra  of  a  PE  oligomer  is  more 
complicated  than  the  spectra  of  a  PE  polymer  chain  because  the  number  of 
degrees  of  freedom  is  much  larger  in  an  oligomer  than  an  infinite,  periodic 
chain.  A  gas  phase  molecule  with  n  atoms  has  3n  -  6  vibrational  modes. 
Thus  the  C6H14  oligomer  used  for  comparisons  with  ab-initio  calculations 


► 

has  54  (3x-6)  vibrational  modes  in  its  spectra.  The  infinite  PE  single  chain 
I  has  14  (3x6-4)  genuine  vibrations  corresponding  to  the  six  atoms  in  the 

C2H4  repeat  unit,  where  three  translations  and  one  rotation  about  the  chain 

axis  do  not  result  in  vibrations.  14  Many  of  the  vibrational  modes  in  the 

*  oligomer  are  nearly  the  same,  but  because  the  oligomer  is  finite,  the 
energies  of  these  modes  are  not  degenerate  as  they  are  in  the  infinite 
polymer.  The  experimental  and  calculated  vibrational  frequencies  for  the 
C6H14  oligomer  are  listed  in  Table  3.  Both  the  observed  frequencies  and 

valence  force  field  (molecular  mechanics)  calculated  frequencies  are  taken 
I  from  Snyder  and  Schachtschneider  .28, 29  The  behavior  of  the  vibrational 

spectra  of  the  n-CH2  series  is  also  thoroughly  discussed  in  these  references. 
It  is  seen  that  the  valence  force  field  calculation  agrees  very  well  with  the 

•  observed  spectrum,  which  is  not  too  surprising  since  this  approach  uses 
force  constants  which  are  obtained  from  spectra.  The  SE  calculated 
frequencies  tend  to  be  systematically  higher  than  observed.  The  ab-initio 
calculated  frequencies  agree  better  with  the  observed  frequencies;  however, 
the  agreement  is  not  as  good  as  can  be  obtained  using  ab-initio  methods 

because  electron  correlation  was  not  included  in  this  calculation.  Even 

i 

without  correlation,  this  ab-initio  calculation  required  in  excess  of  250 
hours  on  an  Elxsi  mini-supercomputer,  while  the  SE  calculation  only  took 
^  about  16  minutes  on  a  Silicon  Graphics  workstation. 

I 

» 
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TABLE  3 


COMPARISON  OF  OBSERVED  AND  CALCULATED  C6H14 
VIBRATIONAL  FREQUENCIES 


OBSERVED 

[cm  -i] 

CALCULATED  [cm-i] 

Valence  Force  Field 

SE  AMI 

Ab-Initio 

— 

61 

55 

-244 

— 

94 

71 

-243 

— 

125 

101 

81 

— 

139 

156 

107 

— 

208 

160 

135 

— 

216 

166 

170 

— 

303 

336 

312 

373 

370 

415 

374 

— 

474 

509 

495 

721 

723 

747 

746 

— 

740 

788 

786 

798 

798 

863 

895 

886 

887 

954 

922 

— 

894 

989 

977 

896 

896 

1013 

993 

996 

1000 

1028 

1027 

1010 

1009 

1182 

1062 

1041 

1036 

1197 

1070 

1067 

1064 

1198 

1103 

1060 

1065 

1207 

1137 

1143 

1145 

1219 

1190 

— 

1178 

1225 

1245 

1225 

1226 

1226 

1273 

1242 

1242 

1237 

1318 

— 

1277 

1261 

1346 
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TABLE  3  (continued) 


OBSERVED  [cm  -i]  CALCULATED  [cm-i] 

Valence  Force  Field  SE  AMI  Ab-Initio 


1302 

1302 

1264 

1371 

1302 

1303 

1312 

1378 

1302 

1303 

1369 

1391 

1353 

1356 

1394 

1431 

— 

1368 

1394 

1439 

1370 

1372 

1396 

1495 

— 

1397 

1396 

1496 

1450 

1444 

1405 

1559 

1452 

1448 

1408 

1564 

1462 

1458 

1408 

1571 

1463 

1462 

1411 

1573 

— 

1462 

1432 

1574 

1463 

1463 

1439 

1577 

— 

1470 

1444 

1586 

1475 

1474 

1458 

1589 

2850 

2852 

3002 

3024 

2851 

2855 

3012 

3030 

— 

2858 

3025 

3048 

2871 

2862 

3037 

3050 

2875 

2882 

3062 

3059 

2885 

2882 

3062 

3061 

— 

2915 

3064 

3062 

2907 

2919 

3064 

3073 

2920 

2924 

3078 

3093 

2934 

2929 

3088 

3103 

2965 

2965 

3098 

3137 

2965 

2965 

3107 

3137 

2965 

2965 

3157 

3148 

2965 

2965 

3157 

3149 
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Using  the  cluster  method,  the  spectra  is  calculated  as  though  the  cluster 
represented  one  repeat  unit  of  a  polymer.  Thus,  for  PE,  the  cluster  method 
assumes  that  the  repeat  unit  is  C6H12  and  not  the  correct  C2H4.  Thus  51, 

3n-3,  vibrational  frequencies  are  calculated  (the  mode  corresponding  to 
rotation  around  the  chain  axis  is  not  discarded  in  the  SE  calculation).  If  a 
one  repeat  unit  cluster  were  long  enough  to  meet  the  criterion  for 
noninteraction  for  the  cluster  method  discussed  earlier,  the  calculated 
spectra  could  be  more  readily  matched  with  the  measured  chain  spectra. 

Since  the  periodic  boundary  conditions  require  that  the  PE  cluster 
contain  at  least  3  repeat  units,  the  vibrational  spectrum  for  the  cluster  is 
much  more  complicated  than  the  PE  polymer  spectrum.  Many  of  the 
cluster  vibrations  involve  C-H  bends  and  stretches,  which  should  not  be 
very  sensitive  to  changes  in  the  backbone  structure  of  PE  caused  by  strain, 
and  will  have  a  minor  effect  on  the  chain  modulus.  In  order  to  make  the 
analysis  tractable  (for  this  study),  the  strain  dependent  shift  of  only  those 
vibrational  modes  which  are  directly  involved  in  backbone  deformation 
(C-C  stretch  modes)  and  can  be  unambiguously  matched  with  observed  PE 
chain  vibrational  modes  (by  matching  normal  coordinates  )  will  be 
discussed  in  detail. 

In  the  MOPAC  calculation,  one  specifies  the  length  of  the  (PE)  cluster, 
thus  one  calculates  the  strain  dependent  properties  of  the  cluster,  here  the 
strain  dependent  frequency  shifts.  Experimental  frequency  shifts  are 
reported  either  as  a  function  of  strain  or  stress,  but  typically  do  not  report 
stress-strain  curves.  The  MOPAC  calculated  strain  dependent  frequency 


shifts  can  be  converted  to  stress  dependent  shifts  using  the  calculated  force 
curve  (first  derivative  of  the  heat  of  formation  potential)  to  determine  the 
proportionality  between  calculated  strain  and  stress. 

Figure  12a  shows  the  stress  dependent  fr  equency  shift  of  the  two 
symmetric  Raman  active  C-C  stretch  modes  in  the  3  repeat  unit  PE  cluster. 
(In  the  infinite  chain  these  modes  are  degenerate.)  The  observed 
equilibrium  (zero  strain  /  zero  stress)  frequency  for  this  mode  is  1059  cm-1, 
while  the  SE  calculated  frequency  is  1157  cm-1.  The  calculated  shift  of 
these  lines  is  linear  out  to  a  stress  of  14  GPa  (4%  strain)  with  a  slope  of -5.86 
cm1  /  GPa  (—  23.3  cm-1  /  %  strain).  Fig  12b,  reproduced  from  reference  15, 
shows  the  measured  stress  dependent  shift  of  this  C-C  stretch  mode  in 
ultra-oriented  PE  film.  Wool19  states  that  measured  shift  in  the  low  stress 
region  is  considered  more  representative  of  the  value  obtained  under  the 
assumption  of  uniform  stress  on  the  chains.  Deviation  from  linearity  in  the 
calculated  stress  dependent  shift  is  not  observed  because  stress  is 
(inherently)  uniformly  applied  in  the  calculation  independent  of  the 
amount  of  stress  applied.  Because  strain  is  proportional  to  stress,  the  fact 
that  the  calculated  stress  dependent  frequency  shift  is  0.52  times  the 
measured  shift  means  that  the  calculated  strain  dependent  frequency  shift 
is  1.92  (1  /  0.52)  larger  than  the  strain  dependence  of  the  measured 
frequency  shift. 

Figure  13a  shows  the  calculated  stress  dependent  frequency  shift  of  the 
asymmetric  C-C  stretch  modes  in  PE.  Because  this  is  an  asymmetric 
stretch,  only  one  frequency  is  present  with  the  correct  normal  coordinates. 
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Rgure  12.  PE  asymmetric  C-C  stretch  mode  stress  dependent  shift  (a)  AMI 
calculated  (b)  experimental.  The  calculated  stress  dependent  shift  is  roughly 
one-half  the  experimental  shift,  which  means  the  calculated  strain  dependent  shift 
is  roughly  twice  the  experimental 
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Figure  13.  PE  symmetric  C-C  stretch  mode  stress  dependent  shift  (a)  AMI 
calculated  (b)  experimental.  The  calculated  stress  dependent  shift  is  roughly 
one-half  the  experimental  shift,  which  means  ihe  calculated  strain  dependent  shift 
is  roughly  twice  the  experimental 
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The  experimental  zero  strain  frequency  for  this  mode  is  1127  cm-1,  while 
the  SE  calculated  frequency  is  1224  cm-1  The  calculated  shift  of  this  mode  is 
-  3.3  cm-1  /%  strain,  -2.88  cm-1  /  GPa.  Figure  13b, reproduced  from 
reference  14  shows  a  measured  shift  of  -5.9  cm-1  /GPa.  The  calculated 
strain  dependent  shift  rate  is  2.04  times  the  observed  strain  dependent  shift 
rate. 

A  possible  explanation  for  the  difference  in  the  calculated  and  observed 
strain  dependent  frequency  shifts  might  be  that  the  measured  bulk  strain  is 
not  be  completely  taken  up  in  c  -axis  chain  strain.  If  the  stress  applied  to 
PE  were  roughly  equally  partitioned  between  chain  slippage  and  chain 
extension,  only  one-half  of  the  bulk  strain  would  be  c-axis  strain,  and  the 
calculated  frequency  shifts  would  agree  with  the  measured  shifts.  This 
hypothesis  was  investigated  by  calculating  the  interaction  energy  of  two  PE 
chains  (inter-chain  friction)  packed  in  the  crystal  structure.  The  program 
Cerius  30  for  Silicon  Graphics  workstations  was  used  for  this  crystal 
simulation.  A  standard  Lennard-Jones  potential  with  an  interaction 
radius  of  5A  was  used  to  describe  the  non-bonded  interactions  between  two 
PE  chains.  The  equilibrium  geometry  obtained  from  the  AMI  calculations 
was  used  for  the  PE  chains.  One  PE  chain  was  moved  along  the  c-axis  of 
the  crystal  with  respect  to  the  other,  and  the  change  in  energy  monitored. 
As  expected,  the  interaction  energy  increased  as  the  chains  were  slipped 
against  each  other;  the  energy  was  periodic  with  a  period  equal  to  the 
repeat  unit  length.  The  energy  barrier  to  chain  slippage  was  found  to  be 
approximately  2%  of  the  chain  axis  barrier.  This  means  that  under  strain 


amorphous  PE  will  initially  deform  by  chain  slippage  rather  than  by 
extension  of  the  chains  themselves  .  (  This  behavior  has  been  observed  in 
unoriented  PE.  31)  Because  the  resistance  to  chain  slippage  is  so  much 
smaller  than  the  chain  resistance,  equipartition  between  inter  and  intra 
molecular  potentials  cannot  be  invoked  to  explain  the  discrepancy  between 
the  measured  and  calculated  frequency  shifts. 

Since  the  discrepancy  in  frequency  shifts  cannot  be  explained  by  inter- 
molecular  interactions,  the  next  logical  step  is  to  investigate  the  validity  of 
the  intra-molecular  interactions,  i.e.,  the  bonding  potentials.  Wool  17 
performed  a  analysis  of  the  strain  dependent  frequency  shifts  of  PE.  In 
addition  to  very  thorough  measurements,  he  also  calculated  PE 
deformation  using  a  molecular  dynamics  model.  The  Morse  potential 
shown  in  Figure  14  was  used  to  describe  the  C-C  bonds  in  PE.  In  this 
reference,  the  source  for  the  parameters  in  the  C-C  Morse  potential  is  not 
listed,  however,  these  parameters  are  essentially  identical  to  those  used  in 
similar  models  found  in  other  references24.  Because  Wool's  model 
matches  both  the  equilibrium  frequencies  and  the  strain  dependent 
frequency  shifts  for  all  the  PE  chain  vibrations  which  were  intense  enough 
to  be  reliably  measured,  this  description  of  the  C-C  bond  in  PE  is  judged  to 
be  accurate. 

The  C-C  bond  strain  potential  which  arises  from  the  AMI  parameters 
was  obtained  by  varying  the  C1-C2  bond  distance  in  the  PE  cluster  shown  in 
Fig  5  and  allowing  all  the  other  variables  to  optimize.  This  potential  is  also 
shown  in  Fig  14  from  -25  %  to  50%  strain.  In  this  figure  it  is  seen  that  the 
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Figure  14.  Comparison  of  Morse  and  AMI  C-C  strain  potentials. 
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Figure  15.  Comparison  of  Morse  and  AMI  C-C  strain  potentials: 
expanded  view  about  equilibrium. 
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AMI  C-C  potential  is  systematically  higher  in  energy  than  the  Morse 
potential.  An  expanded  view  around  equilibrium  is  shown  in  Fig  15,  which 
shows  that  there  is  a  significant  difference  in  the  two  potentials  even  for 
small  variations  about  equilibrium.  The  second  derivative  of  the  AMI  C-C 
potential  at  equilibrium  is  580  ±  20  Nt  /  m,  which  is  31%  stiffer  than  the 
second  derivative  of  the  Morse  C-C  potential  at  equilibrium,  440  Nt  /  m. 

The  effective  CCC  bond  angle  potential  can  be  obtained  by  varying  angle 
Cl— C2-C3  (Fig  5)  and  allowing  the  rest  of  the  cluster  to  optimize.  A  value 
of  1.566e-18  Nt  m  /  rad2  was  obtained  for  the  angle  bending  force  constant 
calculated  by  the  AMI  Hamiltonian.  A  force  constant  for  a  CCC  angle  bend 
of  1.084e-18  Nt  m/rad2  was  used  in  reference  27.  Again,  the  value  from 
AMI  is  approximately  30%  higher  than  has  been  empirically  found  to 
predict  the  correct  frequencies  and  frequency  shifts. 

Wool  reports  a  value  of  267  GPa  for  the  chain  modulus  of  PE  based  on  his 
measurements  and  modelling,  which  differs  by  -10%  from  the  "accepted" 
value  of  300  GPa  for  PE.  The  AMI  result  was  400  GPa,  which  differs  from 
the  "accepted"  value  by  +33%.  The  frequency  shifts  predicted  by  AMI  are 
roughly  a  factor  of  two  larger  than  experimentally  observed.  The  Morse 
potential  used  by  Wool  has  second  and  third  derivatives  that  agree  witn 
experiment.  The  C-C  bond  potential  generated  by  the  AMI 
parameterization  is  judged  to  be  inaccurate  because  its  derivatives  a  re 
larger  than  experimentally  observed,  hence  the  modulus,  which  depends 
directly  in  the  second  derivatives  of  the  potentials,  is  judged  to  be 
inaccurate;  however  it  is  premature  1  o  assign  a  quantitative  error  in  the 
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V  POLYDIACETYLENE  RESULTS  AND  DISCUSSION 
Substituent  Effects 

Polydiacetylene  (PDA)  is  one  of  a  few  polymers  which  can  be 
polymerized  as  a  single  crystal.  The  backbone  chain  structure,  shown  in 
figure  16,  is  simple;  however,  the  substituents  (indicated  by  R  in  figure  16), 
are  typically  large  structures  containing  up  to  30  atoms  (including 
hydrogens).  If  one  considers  a  "ball  and  spring"  mechanical  model  of 
PDA,  it  seems  reasonable  to  assert  that  the  stiffness  of  the  polymer  along  its 
c-axis  should  be  independent  of  substituent  because  the  substituent  spring 
(bond  strength)  is  orthogonal  to  the  chain  axis.  In  reality,  the  bond  between 
the  backbone  and  the  substituent  will  affect  the  bonding  on  the  backbone, 
but  unless  the  backbone  is  bonded  to  a  halogen  atom,  the  influence  on 
backbone  bonding  should  be  small.  One  could  argue  that  the  conjugated 
PDA  backbone  determines  the  chain  stiffness;  however,  because  the 
substituent's  size  determines  the  crystal  packing  density  hence  the  unit  cell 
cross  sectional  area,  substituent  geometry  ultimately  determines  the 
modulus. 

If  this  argument  were  valid,  the  PDA  derivative  which  has  the 
substituents  which  give  the  best  crystal  packing  (highest  density)  would 
have  the  largest  modulus.  Table  4  shows  a  comparison  of  4  PDA  variants 
whose  chemical  formulae  are  listed  in  figure  16.  (PDA  derivatives  are 
commonly  referred  to  by  only  the  substituent  type;  for  clarity  in  this  thesis  I 
will  add  the  abbreviation  PDA.  For  example,  PDA-EUHD  referred  to 
herein  is  referred  to  simply  as  EUHD  in  the  literature).  PDA-EUHD, 


R  R 


PUHD  R  =  CH  OCONHC  H 

2  6  5 

TSHD  R  =  CH2-0-S02C6H4CH3 

O 

II 

EUHD  R  =  CH-0-C-N-C2H5 

H 

DCHD  R  =  CH-N-(Cg  H4)2 

Figure  16.  PDA  backbone  structure  with  formulae  for 
substituents  discussed  in  the  text. 
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which  has  the  smallest  unit  cell  area  also  has  the  largest  mechanical 
modulus.  Because  all  the  subtituents  listed  have  a  carbon  atom  bonded  to 
the  backbone  their  effect  on  the  backbone  bonding  should  be  the  same.  The 
last  column  in  Table  4  lists  the  product  of  the  modulus  and  the  derivative’s 
cross  sectional  area,  the  chain  restoring  force  which  is  the  chain  stiffness 
times  the  repeat  unit  length.  The  restoring  force  is  seen  to  be  independent 
of  substituent  to  within  experimental  accuracy. 


TABLE  4 

COMPARISON  OF  PDA  DERIVATIVE  PROPERTIES32 


PUHD 

97.4 

45  M 

4.38 

TSHD 

95.5 

43  M 

4.11 

50  R 

4.75 

DCHD 

106.3 

45  M 

4.78 

47  B 

4.95 

EUHD 

65.1 

61 M 

3.97 

74  R 

4.79 

M:  Mechanical  measurement  of  modulus 

R:  Raman  measurement  of  modulus  from  acoustic  velocity  in 

polymer 

B:  Brillioun  measurement  of  modulus 


Ref  32  does  not  give  the  uncertainties  in  the  modulus  measurements  which 
are  taken  from  a  number  of  experiments  (which  are  referenced  in  this 
paper),  but  an  uncertainty  of  ±  5  GPa  is  typically  quoted  in  similar 


experiments. 

The  weak  dependence  of  PDA  modulus  on  substituent  is  very  important 
for  this  study  because  for  every  4  carbon  atoms  in  the  backbone  repeat  unit, 
there  are  as  many  as  50  substituent  atoms.  The  vast  majority  of  the 
computation  time  would  be  spent  optimizing  the  geometry  of  the 
substituents,  and  the  number  of  calculated  vibrational  frequencies  would 
grow  very  large  and  needlessly  complicate  the  analysis  of  the  spectrum 
since  the  substituent  frequencies  would  be  interspersed  with  the  strain 
dependent  backbone  frequencies.  More  importantly,  since  the 
computational  cost  scales  as  n2  (number  of  electrons),  the  MOPAC 
calculations  would  have  been  too  slow  to  be  useful.  Instead  of  calculating 
the  real  polymers,  PDA  with  hydrogens  as  substituents  was  considered  in 
the  majority  of  the  work  performed  in  this  study.  This  is  a  good  model  for 
comparison  because  molecular  mechanics  models  of  PDA  derivatives  often 
treat  the  substituents  as  a  single  atom  with  a  mass  equal  to  the  mass  of  the 
entire  substituent  in  order  to  obtain  a  good  match  of  the  vibrational  coupling 
to  the  backbone. 33.34, 35,36  in  order  to  investigate  isotope  effects,  MOPAC  5.0 
allows  the  user  to  specify  the  masses  of  the  atoms,  e.g.,  specifying  a 
hydrogen  atom  with  mass  2.0  changes  the  hydrogen  to  deuterium.  MOPAC 
5.0  does  not  restrict  the  user  to  naturally  occurring  isotopes;  an  arbitrary 
mass  can  be  assigned  to  any  atom.  This  flexibility  was  exploited  to  match 
the  empirical  models  by  attaching  very  heavy  hydrogen  atoms  to  PDA.  For 
example,  to  match  the  model  of  PDA-HHD  used  in  reference  33,  mass  31.26 
atomic  mass  units  (amu)  and  15.12  hydrogen  atoms  were  used  as 


substituents.  The  mass  of  the  atoms  does  not  affect  the  calculated 
minimum  energy  geometries  because  in  the  Hartree-Fock  method  the 
solution  of  the  electronic  wavefunctions  is  independent  of  the  nuclear 
wavefunctions.  This  was  verified  for  MOPAC  5.0  by  calculating  the 
optimum  geometry  of  a  simple  molecule  while  varying  the  masses  of  the 
hydrogen  atoms  in  the  molecule.  Even  with  mass  106  amu  hydrogen,  all  of 
the  variables,  bond  lengths,  angles,  etc,  were  exactly  the  same  (to  machine 
precision). 

Molecular  Structure 

An  SE  polymer  molecular  structure  calculation  on  PDA  has  not  been 
previously  performed,  so  the  first  step  in  this  study  was  to  determine  the 
minimum  cluster  length.  The  minimum  cluster  length  was  determined  by 
examining  the  heat  of  formation  per  repeat  unit  and  c-axis  length  as  a 
function  of  the  number  of  repeat  units  in  the  cluster.  The  results  are 
summarized  in  Table  5.  These  results  indicate  that  the  normalized  heat  of 
formation  and  c-axis  length  reach  their  asymptotic  values  with  a  4  repeat 
unit  cluster,  which  is  used  for  all  polymer  calculations.  For  comparison, 
Table  6  lists  similar  results  for  the  PDA  oligomer  sequence.  The  oligomers 
demonstrate  behavior  similar  to  the  cluster  sequence.  The  decrease  in  heat 
of  formation  per  repeat  unit  from  1  to  2  occurs  because  an  extra  carbon  is 
included  in  the  one  repeat  emit  case  to  match  the  polymer  c-axis. 


TABLE  5 


PDA  CLUSTER  SERIES 


n:  C4H2 

Heat  of 

Formation  (AHf) 
[kCal  s/mole] 

AHf /n 
[kCal  s/mole] 

Cluster 
Length  [A] 

Repeat  Unit 
Length  [A] 

1 

77.68 

77.68 

4.708 

4.708 

2 

127.21 

63.60 

9.753 

4.876 

3 

186.51 

62.17 

14.634 

4.878 

4 

247.94 

61.98 

19.505 

4.876 

5 

309.83 

61.96 

24.382 

4.876 

6 

371.76 

61.96 

29.282 

4.880 

7 

433.66 

61.95 

34.122 

4.874 

Formula 

TABLE  6 

PDA  OLIGOMER  SERIES 

n  Repeat  Units  AHf  [kCals/mole] 

AHf  (n)-AHf(n-l) 

[kCals/mole] 

c5h6 

1 

56.30 

C8  Hg 

2 

95.92 

39.62 

C12H10 

3 

158.55 

62.63 

C16H12 

4 

220.59 

62.04 

C20H14 

5 

282.53 

61.94 

C24H16 

6 

344.46 

61.93 

C28Hi8 

7 

406.40 

61.94 

A  4  repeat  unit  oligomer  of  PDA  with  hydrogen  substituents  is  too  large 
even  for  a  low  level  ab-initio  calculation.  A  smaller  PDA  oligomer  shown  in 
figure  17  was  used  for  comparisons  between  ab-initio  and  SE  calculations. 
Also  listed  are  geometric  values  from  a  4  repeat  unit  cluster  with  methyl 
(CH3)  substituents  and  an  experimental  X-ray  diffraction  determination  of 
PDA-THD37. 

There  are  several  notable  features  in  the  geometry  of  this  PDA  oligomer. 
First,  the  minimum  energy  geometry  occurs  with  the  end  carbons  both 
above  (cis)  the  plane  of  the  central  acetylene  group,  while  in  the  polymer 
successive  repeat  units  are  always  on  opposite  sides  of  this  plane  (t.rans). 

In  the  oligomer  with  hydrogen  substituents,  the  barrier  to  rotation  was 
calculated  and  found  to  be  very  low  (thermal),  which  means  that  the  ends  of 
this  oligomer  would  freely  rotate  at  room  temperature.  Steric  interactions 
with  larger  substituents  would  restrict  the  rotation  of  the  oligomer; 
however, in  the  AMI  cluster  calculations  with  hydrogen  substituents,  the 
PDA  cluster  geometry  always  optimized  in  the  trans  configuration,  which 
is  observed  experimentally  for  the  PDA  polymer.  This  free  rotation  is  a 
significant  end  effect  in  the  oligomer,  and  complicated  attempts  to  calculate 
the  modulus  on  a  two  repeat  unit  oligomer.  (The  rationale  for  these 
attempts  is  discussed  later  in  this  chapter.) 

There  are  significant  differences  in  the  backbone  C-C  bond  lengths 
between  PDA  and  PE  .  The  PDA  C-C  single  bond  is  considerably  shorter 
than  both  the  C-C  single  bond  in  PDA,  and  the  C-C  single  bond  to  the  carbon 
atom  in  the  substituent.  This  occurs  because  the  backbone  "single"  C-C 
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Figure  17.  PDA  geometry,  calculated  for  the  oligomer  shown,  for  a  4  repeat  unit  cluster  with 
hydrogen  and  methyl  substituents,  and  values  obtained  from  X-ray  diffraction  measurements. 


bond  is  really  a  mixture  of  a  single  C-C  bond  and  a  triple  C-C  bond.  This 
bond  is  somewhere  between  a  pure  C-C  single  bond  with  a  length  of  1.53  A 
and  a  pure  C-C  double  bond  with  a  length  of  1.326  A.  The  AMI  calculation 
predicts  that  the  C-C  double  bond  (C1-C2)  and  C-C  single  bond  (C2-C3)  in 
PDA  have  nearly  equal  lengths;  the  ab-initio  calculation  shows  these  bond 
lengths  to  differ  measurably,  which  agrees  much  better  with  experiment 
than  the  SE  results.  The  error  in  these  AMI  bond  lengths  is  larger  than  the 
average  error  in  bond  lengths  as  listed  in  Table  2,  which  may  show  a 
limitation  in  the  AMI  parameterization  for  conjugated  carbon  bonds.  The 
C— C  triple  bond  lengths  calculated  by  AMI  and  ab-initio  both  agree  very 
closely  with  experiment.  The  bond  angles  calculated  by  AMI  are  within  the 
published  accuracy  of  this  method.  An  earlier  MNDO  calculation  on  a  one 
repeat  unit  PDA  oligomer  with  methyl  substituents  showed  results  which 
were  more  consistent  with  the  experimentally  observed  geometry.38  A 
MNDO  calculation  on  the  oligomer  shown  in  Fig  17  was  performed  in  this 
study  and  found  to  produce  approximately  the  same  results  as  the  AMI 
calculation. 

Figure  18  shows  a  comparison  of  the  distribution  of  strain  within  the 
PDA  repeat  unit  as  a  function  of  c-axis  strain.  The  bond  angle  is  seen  to 
show  the  largest  relative  change  with  strain.  Initially  the  single  and 
double  C-C  bonds  deform  at  the  same  relative  rate,  but  at  strains  greater 
than  4%  the  single  C-C  bond  begins  to  change  more  rapidly.  The  triple  C-C 
deforms  very  little  at  low  strains  and  begins  to  deform  linearly  for  strains 
above  2%.  Overall,  this  behavior  is  consistent  with  what  one  would  expect 
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Figure  18.  PDA  distribution  of  chain  strain  among  bond  angle 
and  bond  lengths. 
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based  on  the  anticipated  stiffness  of  ea^h  mode,  the  triple  bond  being  the 
strongest,  the  angular  resistance  being  the  weakest;  however,  it  is 
interesting  that  the  relative  deformation  of  the  PDA  single  and  double  C-C 
bonds  is  roughly  equal.  This  result  may  be  explained  by  the  fact  that  these 
are  not  pure  double  and  single  bonds.  The  AMI  calculated  deformation  is 
consistent  with  the  experimentally  reported  deformation  of  PDA-HDU, 
(R=CH2OCONHC6H5).39 

Modulus  Calculations 

The  strain  dependent  heat  of  formation  curves  for  a  4  repeat  unit  PDA 
cluster  is  shown  in  figure  19  from  -5%  to  10  %  strain.  Two  cases  are  shown 
in  this  figure,  one  with  hydrogen  substituents  and  one  with  methyl 
substituents.  As  is  seen,  these  potentials  are  essentially  identical,  which 
reinforces  the  argument  that  PDA  chain  stiffness  is  largely  independent  of 
substituent.  Results  below  -5%  were  not  obtained  because  the  specification 
of  the  PDA  geometry  caused  the  optimization  to  fail  if  a  bond  torsion  angle 
rotation  occurred.  Numerous  alternate  geometry  specifications  were 
attempted  to  work  around  this  problem,  but  none  were  successful.  Since 
experimental  data  below  -4%  have  not  been  reported,  this  issue  was  not 
pursued.  The  chain  stiffness  obtained  by  the  same  procedure  used  for  PE  is 
42  Nt/m.  The  AMI  calculated  modulus  for  the  PDA  variants  described 
before  are  listed  in  Table  7  along  with  the  measured  bulk  values  (cross 
sectional  area  are  taken  from  Table  4).  As  was  seen  in  PE,  the  AMI 
calculated  modulus  is  substantially  higher  than  the  measured  values. 
Better  agreement  was  anticipated  because  PDA’s  structure  is  much  closer 
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Figure  19.  Comparison  of  Mopac  heat  of  formation  curve  for  a  4  repeat 
unit  PDA  cluster  with  hydrogen  and  methyl  substituents 
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to  the  ideal  chain  simulated  in  these  calculations.  As  was  required  for  PE, 
other  predictions  of  the  AMI  calculation  need  to  be  compared  with 
experiment  to  assess  the  accuracy  of  this  method. 


TABLE  7 

COMPARISON  OF  EXPERIMENTAL  AND  CALCULATED  PDA  CHAIN 

MODULUS 

PDA  Experimental  Modulus  |Gpa]  32  AMI  Modulus  [GPa] 

Derivative 


PUHD 

45  M 

86 

TSHD 

43  M 

50  R 

88 

DCHD 

45  M 

47  B 

79 

EUHD 

61 M 

74  R 

129 

In  order  to  discover  how  much  of  the  modulus  error  comes  from  the 
approximations  made  in  the  AMI  method,  a  direct  comparison  of  the 
strained  geometry  heat  of  formation  calculations  between  the  SE  and  ab- 
initio  methods  was  planned.  The  ab-initio  codes  which  was  available  for 
this  study,  Gaussian88  and  Gaussian90,  do  not  have  provision  for  the  use  of 
periodic  boundary  conditions  to  simulate  an  infinite  chain,  so  this 
comparison  had  to  be  limited  to  an  oligomer.  Because  of  the  memory 
constraints,  the  maximum  practical  oligomer  which  could  be  used  in  the 
ab-initio  calculation  was  only  two  repeat  unit  long.  This  oligomer  is  shown 
in  figure  20.  Calculations  were  first  performed  using  the  SE  method  to  see 
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if  the  effect  of  strain  on  a  two  repeat  unit  oligomer  was  representative  of  the 
infinite  chain.  Two  interesting  results  were  obtained  from  the  SE 
calculations.  First,  the  vector  direction  in  which  strain  is  applied  affects 
the  optimized  geometry.  Strain  is  applied  in  an  oligomer  by  displacing  the 
"end"  atoms  (symmetrically)  a  small  distance  away  from  their  equilibrium 
location  and  holding  them  fixed.  The  rest  of  the  oligomer  geometry  is 
allowed  to  optimize  under  this  constraint.  Because  the  atom  locations  are 
specified  by  a  bond  distance  and  angle,  successive  displacements  take  place 
along  an  effective  strain  vector.  Several  vector  directions  which  were  used 
are  shown  in  figure  20.  Because  each  successive  atom  displacement  is 
very  small,  the  actual  position  where  the  atom  is  placed  differs  minutely 
depending  on  the  strain  vector.  Surprisingly  (at  least  to  the  author),  the 
heat  of  formation  and  optimized  geometry  was  found  to  be  very  sensitive  to 
the  strain  vector.  In  order  to  match  the  strain  direction  in  the  cluster 
method,  atom  displacement  was  specified  along  vector  III,  which  is  along 
the  polymer  c-axis.  Unfortunately,  the  results  of  this  calculation  showed 
that  a  two  repeat  unit  oligomer  was  dominated  by  end  effects,  and  the 
comparison  to  t^e  infinite  chain  was  invalid.  Therefore,  the  ab-initio 
comparison  was  not  performed. 

Comparisons  With  Spectroscopy 

The  accuracy  of  the  modulus  calculation  of  PDA,  as  for  PE,  can  be 
investigated  by  comparing  the  measured  and  AMI  calculated  strain 
dependent  vibrational  frequency  shifts.  The  large,  4  repeat  unit  cluster 
required  for  PDA  made  this  analysis  more  difficult  because  the  normal 


vibrational  modes  for  a  4  repeat  unit  cluster  are  very  different  than  the 
vibrational  modes  for  the  infinite  chain.  Because  the  substituents  couple 
into  some  of  the  backbone  chain  Raman  active  modes,  it  is  difficult  to 
accurately  match  these  modes  with  the  hydrogen  substituents.  The 
experimental  spectrum  is  complicated  because  some  frequencies  that  arise 
solely  from  the  substituents  are  interspersed  with  the  chain  modes.  The 
substituent  frequencies  do  not  shift  with  strain,  so  they  don't  provide  any 
information  on  the  accuracy  of  the  strain  dependence  of  the  AMI 
calculation.  After  a  thorough  analysis,  only  two  calculated  vibrational 
frequencies  could  be  unambiguously,  reliably  matched  with  the  reported 
spectra.  The  normal  coordinates  of  these  modes  are  shown  in  figure  21. 
34,36  These  measurements  have  been  performed  on  a  variety  of  PDA 
substituents  of  varying  degrees  of  polymerization,  so  there  is  considerable 
spread  in  the  reported  frequency  shifts.40 

The  strain  dependent  frequency  calculations  were  first  performed  using 
mass  1  amu  hydrogen  substituents,  which  produced  an  unrealistic 
spectrum  because  the  coupling  of  the  light  mass  hydrogens  to  the  backbone 
vibrations  was  too  strong.  When  the  calculation  was  performed  with  mass 
32  hydrogens,  the  backbone  coupling  better  matched  the  experiment  and 
the  spectrum  was  easier  to  interpret.  Interestingly,  the  frequencies  of  the 
two  strongly  strain  dependent  Raman  modes  were  found  to  be  independent 
of  the  hydrogen  mass.  A  possible  problem  with  approximating  the 
substituents  with  heavy  hydrogen  atoms  is  that  the  bond  length  and  the 
force  constant  of  the  atom  bonded  to  the  backbone  are  wrong,  which  might 
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Figure  21.  Eigenvectors  for  the  two  PDA  Raman  active 
backbone  vibrational  modes  whose  frequency  shifts  are 
investigated  in  this  study 
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affect  the  vibrational  coupling  to  the  backbone.  This  possibility  was 
investigated  by  performing  a  frequency  calculation  on  PDA  with  methyl 
substituents,  which  gives  the  C-C  bond  to  the  backbone  found  in  all  the  PDA 
derivatives,  but  changing  the  mass  of  the  carbon  atom  on  the  methyl  group 
to  32  amu.  The  backbone  frequencies  described  before  did  not  change.  This 
was  a  useful  result  for  this  study  because  the  calculation  time  for  a  methyl 
substituted  PDA  cluster  is  very  much  longer  than  hydrogen  substituted 
PDA  because  the  methyl  substituted  cluster  has  27  more  (hydrogen)  atoms. 

The  strain  dependent  frequency  shift  of  the  backbone  frequency 
corresponding  the  mode  shown  in  figure  21a  is  shown  in  figure  22. 

Because  a  four  repeat  unit  cluster  was  used  in  the  calculation,  there  are 
four  frequencies  in  the  cluster  which  correspond  to  the  vibrational  mode  in 
the  infinite  chain.  The  difference  in  these  frequencies  is  very  small,  which 
shows  that  they  are  nearly  degenerate.  The  calculated  frequency  shift  as  a 
function  of  strain  shows  more  spread  than  the  PE  results,  which  is  a 
limitation  of  the  accuracy  in  calculating  the  frequencies  on  such  a  large 
cluster.  The  calculated  strain  dependent  frequency  shift  of  this  mode  is 
40±3  cm-i  /  %  strain,  which  is  the  average  of  the  four  lines.  The 
experimental  value  of  the  frequency  shift  varies;  but  the  largest  reported  is 
20  ±  0.5  cm-i  /%  strain.32  The  experiment  that  measured  the  largest  shift 
was  performed  on  PDA  samples  with  best  morphologies,  which  implies 
that  this  shift  is  more  representative  of  the  defect  free  crystal.  Galliotis  et 
al2i  list  a  number  of  measured  values  for  the  shift  of  two  backbone  Raman 
active  modes  for  a  range  of  PDA  derivatives  of  varying  quality,  which  gives 
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Figure  22.  Strain  depenedent  frequency  shifts  of  the  4  repeat  unit  PDA 
cluster  corresponding  to  the  vibrational  mode  shown  in  figure  20  (a). 
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an  idea  of  the  sensitivity  of  these  shifts  to  sample  preparation.  As  was 
found  for  PE,  the  calculated  frequency  shift  of  this  mode  is  approximately 
twice  that  observed  experimentally. 

The  frequency  shift  of  the  mode  shown  in  figure  21b  is  shown  in  figure 
23.  There  are  only  3  frequencies  in  the  4  repeat  unit  cluster  for  this  mode 
because  it  involves  displacement  of  the  carbons  along  the  double  bonds,  and 
there  are  only  3  double  bonds  in  the  cluster.  There  is  considerable  spread  in 
these  calculated  strain  dependent  frequencies;  however,  the  average 
calculated  frequency  shift  of  this  mode  is  15  cm-1  /  %  strain  ,  which  is 
approximately  twice  the  highest  reported  value  for  the  shift  of  this  mode15. 

As  was  the  case  for  PE,  it  would  be  informative  to  compare  the  potentials 
generated  by  AMI  for  the  individual  bonds  in  PDA  with  Morse  potentials; 
however,  Morse  potential  parameters  have  not  been  reported  for  these 
bonds.  Because  the  electrons  are  not  localized  to  specific  atom  pairs  in 
PDA  ,it  is  more  difficult  to  define  potentials  that  accurately  describe  the 
interaction.  The  AMI  potentials  for  the  single,  double,  and  triple  C-C  bonds 
in  PDA  are  shown  in  Figure  24.  These  potentials  are  best  described  by 
third  order  polynomials,  which  are  stiffer  than  Morse  potentials.  Instead 
of  explicitly  defining  a  bond  potentials,  empirical  force  constant  approaches 
for  calculating  the  modulus  in  various  PDA  derivatives  obtain  set  of  force 
constants  from  the  equilibrium  spectrum,  and  use  a  fudge  factor  to  describe 
the  variation  of  the  force  constant  with  strain.  An  empirically  derived 
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Figure  23.  Strain  depenedent  frequency  shifts  of  the  4  repeat  unit  PDA 
cluster  corresponding  to  the  vibrational  mode  shown  in  figure  20  (b). 
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Single  B< 


estimate  of  the  force  constant  variation  is  33>34 


r-A  (14) 

where  A=1.06  A  and  B  =  1.65  x  10-®  Nt.  This  approximate  relation  is  known 
to  overpredict  the  strain  dependent  frequency  shift,  which  is  not  surprising 
because  it  predicts  the  force  constant  to  fall  off  as  1/r,  which  is  slower  than 
the  force  constant  variation  of  a  Morse  potential.  An  accurate  calculation  of 
the  strain  dependent  frequency  shift  in  PDA  has  not  been  reported. 

A  comparison  of  the  equilibrium  force  constants  obtained  from  the 
second  derivatives  of  the  AMI  potentials  shown  in  Fig  24  and  those  reported 
in  by  Batchelder  and  Wu  is  shown  in  Table  8.  As  is  seen,  the  AMI  values 
for  the  force  constants  are  considerably  higher  than  those  derived  from  the 
PDA  spectra.  The  models  from  which  these  PDA  experimental  force 
constants  were  derived  were  simplified,  so  these  force  constants  are  not 
guaranteed  to  be  unique.  Also  listed  in  Table  8  are  the  force  constants  for 
the  corresponding  unconjugated  C-C  bonds  calculated  by  AMI.  The  force 
constants  for  the  unconjugated  bonds  are  in  a  ratio  which  is  close  to  1:2:3 
for  single  to  double  to  triple,  while  in  PDA  the  single  and  double  have 
roughly  equal  strength,  each  approximately  one-half  that  of  the  triple  bond. 


TABLE  8 


COMPARISON  OF  PDA  BACKBONE  FORCE  CONSTANTS  [Nt/m] 


C-C  Bond 
Type 

AMI 

Ref  29 

Ref  30 

AMI 

Unconjugated 

Single 

820 

450 

560 

620  (C2H6) 

Double 

880 

550 

510 

1220  (C2H4) 

Triple 

1530 

1200 

1170 

1880  (C2H2) 

The  deformation  of  PDA  does  not  involve  only  one  kind  of  bond,  so  a 
comparison  of  the  individual  calculated  and  experimental  force  constants 
does  not  have  a  direct  correlation  to  the  error  of  the  calculated  modulus. 
However,  as  a  rough  comparison,  using  a  simple  addition  of  the  bonds 
found  in  the  repeat  unit,  two  single  C-C,  one  double  C-C,  and  one  triple  C-C, 
bond,  the  AMI  calculated  sum  is  approximately  40%  higher  than 
experiment.  Because  the  deformation  analysis  showed  that  the  change  of 
the  bond  angle  contributed  more  to  extension  along  the  c-axis  than  bond 
extension,  the  overestimate  of  the  modulus  is  less  than  the  40% 
overestimate  of  the  bond  stiffnesses. 
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VI  SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 

The  ultimate  mechanical  properties  of  polymers,  assuming  perfect 
morphology,  will  be  limited  by  the  mechanical  properties  of  a  single,  ideal 
chain.  The  standard  approach  used  to  predict  the  chain  moduli  of 
polymers,  the  ubiquitous  method  of  Treolar*!,  uses  force  constants  for  bond 
stretching,  bond  angle  opening,  and  torsion  angle  twisting,  obtained  from 
spectroscopy  to  build  up  a  classical  mechanics  model  of  the  polymer  as  a 
series/parallel  combination  of  ideal  springs.  Extensions  to  this  approach 
incorporate  strain  dependence  (anharmonicity)  in  the  force  constants  (ref 
35,  36  for  example).  While  the  method  of  Treolar  has  the  advantage  of  being 
grounded  in  experiment,  it  is  inconvenient  to  apply  in  many  situations 
because  polymer  spectra  tend  to  be  complex,  making  the  task  of 
unambiguously  obtaining  the  required  force  constants  very  difficult.  In 
addition,  the  method  of  Treolar  can  only  be  applied  to  polymers  which  have 
already  been  synthesized,  i.e.,  it  provides  no  predictive  capability  for 
designing  polymers  with  desireable  mechanical  properties.  Likewise, 
molecular  mechanics  is  based  on  situation  specific  experimental  data,  and 
has  a  limited  predictive  capability. 

SE  quantum  mechanical  calculations  of  polymer  mechanical  properties 
were  attempted  because  the  SE  method  is  generic  enough  to  have  a 
predictive  capability,  and  ab-initio  methods  at  present  are  too  time 
consuming  to  be  practical  for  polymers.  Moduli  values  obtained  were 
systematically  (and  significantly)  higher  than  experimental  values,  which 
was  attributed  to  morphological  imperfections  which  will  be  present  in  any 


real  material.  The  goal  of  this  investigation  was  to  explore  possible 
inaccuracies  inherent  in  the  SE  method  which  would  contribute  to  these 
discrepancies.  The  criterion  used  to  make  this  assessment  was  the  strain 
dependent  behavior  of  polymer  spectra  and  strain  potentials  of  individual 
bonds  within  the  polymers  PE  and  PDA,  polymers  which  have  been  well 
characterized  experimentally  and  have  relatively  simple  structures. 

At  the  beginning  of  this  investigation  it  was  anticipated  that  the  SE 
method  would  be  valid  for  some  small  region  around  equilibrium  (whose 
range  would  be  bounded  in  this  investigation)  but  would  be  invalid  above 
(and  below)  some  critical  value  of  strain.  This  anticipation  was  based  on 
the  fact  that  the  SE  method  used  in  the  calculation  (as  well  as  most  SE 
methods)  was  parameterized  to  reproduce  equilibrium  molecular 
properties.  Almost  any  physically  realistic  molecular  potential  is  nearly 
harmonic  around  equilibrium  (harmonic  approximation  of  an  anharmonic 
oscillator)  so  the  extrapolation  of  equilibrium  properties  for  small 
deviations  around  equilibrium  is  valid.  What  was  found  in  this 
investigation  is  that  the  force  constants  obtained  from  the  SE  method,  even 
at  equilibrium,  are  too  high.  The  individual  strain  dependent  bond 
potentials  produced  by  the  SE  method  were  found  to  be  significantly  stiffer 
than  Morse  potentials  used  to  describe  these  bonds;  the  Morse  potentials 
used  in  the  comparison  have  been  shown  to  accurately  describe  the 
anharmonicity  and  reproduce  strain  dependent  spectral  shifts  in  PE.17  As 
a  group,  the  SE  generated  potentials  were  found  to  be  well  fit  by  cubic 
potentials,  which  is  essentially  a  first  order  correction  to  a  harmonic 
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potential.  Because  a  polymer  has  many  deformation  modes,  the 
overprediction  of  bond  stiffness  is  not  necessarily  equal  to  the  overprediction 
of  polymer  modulus. 

The  consequences  of  the  SE  prediction  of  bond  potentials  is  illustrated  in 
Figure  25.  Figure  25a  shows  a  comparison  of  Morse,  harmonic,  and  cubic 
potentials  which  have  been  constructed  to  have  the  same  second  derivative 
at  equilibrium.  The  Morse  potential  is  the  C-C  potential  from  Wool17;  the 
harmonic  and  cubic  potentials  are  least  squares  fits  to  this  Morse  potential. 
Figure  25b  shows  the  first  derivative  of  these  potentials,  which  is  the 
molecular  stress-strain  curve.  Figure  25c  shows  the  second  derivative  of 
these  curves,  which  has  the  same  shape  as  the  modulus  curve  for  this 
potential  because  the  second  derivative  and  modulus  differ  only  by 
constants.  The  second  derivative  of  the  harmonic  curve  is  a  constant;  an 
ideal  harmonic  material  would  have  a  modulus  which  is  independent  of 
strain  (and  would  never  fracture).  If  the  bond  is  accurately  described  by  a 
cubic  potential,  the  modulus  linearly  decreases  with  tension  and  increases 
with  compression.  The  second  derivative  of  the  Morse  potential  initially 
decreases  faster  in  tension  than  the  cubic,  and  falls  off  slower  at  strains 
above  15%.  The  Morse  potential's  second  derivative  increases  faster  than 
the  cubic  in  compression.  However,  for  small  values  around  equilibrium, 
the  second  derivative  of  the  Morse  potential  can  be  well  described  by  a  linear 
fit,  so  it  could  be  difficult  to  discriminate  between  Morse  and  cubic 
potentials  since  experiments  are  limited  to  small  strains  (by  bulk  failure 
mechanisms). 


Figure  25.  Comparison  of  Morse,  Harmonic,  and  Cubic  potentials  (a), 
the  force  curves  (b),  and  the  stiffness  curves  (c)  from  each.  These  potentials 
were  constrained  to  have  equal  second  derivatives  at  equilibrium. 
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What  are  the  possible  reasons  for  the  observed  behavior  of  the  SE 
method?  First,  the  SE  method  was  not  parameterized  with  data  that  give 
information  about  the  derivatives  of  the  molecular  potentials,  even  at 
equilibrium.  Thus  the  inaccuracy  of  the  predicted  stiffness  away  from 
equilibrium  is  not  surprising.  A  more  basic  limitation  may  be  the  fact  that 
SE  methods  are  subject  to  the  Hartee-Fock  limit.  It  is  necessary  to  include 
electron  correlation  in  ab-initio  calculations  to  obtain  good  agreement  with 
experimental  spectra,  i.e.,  derivatives  of  the  molecular  bond  potentials.6 
Electron  correlation  is  known  to  reduce  bond  stiffness;  bonds  calculated  at 
„he  Hartree-Fock  limit  are  found  to  be  systematically  stiffer  than 
experiment,  which  is  consistent  with  the  results  of  the  SE  calculations 
performed  in  this  investigation.  Since  mechanical  properties  of  polymers 
ultimately  depend  on  derivatives  of  the  molecular  bond  potentials,  the 
accuracy  of  a  calculation  of  polymer  chain  modulus  can  be  judged  on  the 
basis  of  the  method's  ability  to  match  these  derivatives  on  systems  where 
experimental  data  exists. 

SE  calculdtions  produce  usefully  accurate  results  for  those  properties 
against  which  they  have  been  benchmarked;  the  AMI  parameterization  did 
not  include  the  spectroscopic  information  which  would  give  better  values 
for  the  equilibrium  force  constant  and  a  more  accurate  description  of  the 
behavior  of  the  system  away  from  equilibrium.  For  the  foreseeable  future, 
ab-initio  methods  will  remain  too  computationally  expensive  to  be  practical 
for  the  prediction  of  the  mechanical  properties  of  polymers.  In  order  to 
improve  the  SE  calculations  of  mechanical  properties  of  polymers,  one  can 
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parameterize  the  method  to  reproduce  spectroscopic  data,  which  will  give  a 
better  anchor  for  bond  stiffness.  Scaling  procedures  to  improve  force 
constants  calculated  at  the  Hartree-Fock  limit  have  been  developed  for 
some  situations.  9  It  may  be  possible  to  incorporate  a  similar  procedure  in 
a  SE  paramterization,  and  avoid  the  computational  cost  of  ab-initio 
methods.  An  extensive  parameterization  could  be  avoided  because  nearly 
all  polymers  of  interest  contain  only  carbon,  hydrogen,  oxygen,  nitrogen, 
sulfur,  and  phosphorous,  which  would  reduce  the  cost  of  such  an  effort. 

The  synthesis  to  high  molecular  weight  (long  chains),  purification,  and 
processing  of  a  new  polymer  is  a  long  and  expensive  process.  The  accurate 
computational  prediction  of  the  ultimate  mechanical  (and  other)  properties 
of  polymers  is  a  desireable  goal  because  it  can  provide  guidance  on  which 
new  polymers  possess  the  best  potential  for  high  modulus  and  strength. 
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